« back to all changes in this revision
Viewing changes to xedit/lisp/mp/mpi.c
- browse files at revision 2
- compare with another revision
- download diff
- download tarball
- view history from revision 2
- Committer: Bazaar Package Importer
- Author(s): Julien Cristau
- Date: 2008-09-23 00:24:45 UTC
- mfrom: (1.1.2 intrepid)
- Revision ID: james.westby@ubuntu.com-20080923002445-mb2rwkif45zz1vlj
Tags: 7.3+4
* Remove xedit from the package, it's unmaintained and broken
(closes: #321434).
* Remove xedit's conffiles on upgrade.
(closes: #321434).
* Remove xedit's conffiles on upgrade.
- files added:
- bitmap
- xgc
- xgc/Written
- files removed:
- debian/xsfbs/xsfbs-autoreconf.mk
- xedit
- xedit/lisp
- xedit/lisp/modules
- xedit/lisp/modules/progmodes
- xedit/lisp/mp
- xedit/lisp/re
- xedit/lisp/test
- files modified:
- debian/changelog
added
removed
xedit/lisp/mp/mpi.c
1
/*
2
* Copyright (c) 2002 by The XFree86 Project, Inc.
3
*
4
* Permission is hereby granted, free of charge, to any person obtaining a
5
* copy of this software and associated documentation files (the "Software"),
6
* to deal in the Software without restriction, including without limitation
7
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
8
* and/or sell copies of the Software, and to permit persons to whom the
9
* Software is furnished to do so, subject to the following conditions:
10
*
11
* The above copyright notice and this permission notice shall be included in
12
* all copies or substantial portions of the Software.
13
*
14
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
17
* THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
18
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
19
* OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
20
* SOFTWARE.
21
*
22
* Except as contained in this notice, the name of the XFree86 Project shall
23
* not be used in advertising or otherwise to promote the sale, use or other
24
* dealings in this Software without prior written authorization from the
25
* XFree86 Project.
26
*
27
* Author: Paulo César Pereira de Andrade
28
*/
29
30
/* $XFree86: xc/programs/xedit/lisp/mp/mpi.c,v 1.12 2002/11/20 07:44:43 paulo Exp $ */
31
32
#include "mp.h"
33
34
#ifdef __UNIXOS2__
35
# define finite(x) isfinite(x)
36
#endif
37
38
/*
39
* Prototypes
40
*/
41
/* do the hard work of mpi_add and mpi_sub */
42
static void mpi_addsub(mpi *rop, mpi *op1, mpi *op2, int sub);
43
44
/* logical functions implementation */
45
static INLINE BNS mpi_logic(BNS op1, BNS op2, BNS op);
46
static void mpi_log(mpi *rop, mpi *op1, mpi *op2, BNS op);
47
48
/* internal mpi_seti, whithout memory allocation */
49
static void _mpi_seti(mpi *rop, long si);
50
51
/*
52
* Initialization
53
*/
54
static BNS onedig[1] = { 1 };
55
static mpi mpone = { 1, 1, 0, (BNS*)&onedig };
56
57
/*
58
* Implementation
59
*/
60
void
61
mpi_init(mpi *op)
62
{
63
op->sign = 0;
64
op->size = op->alloc = 1;
65
op->digs = mp_malloc(sizeof(BNS));
66
op->digs[0] = 0;
67
}
68
69
void
70
mpi_clear(mpi *op)
71
{
72
op->sign = 0;
73
op->size = op->alloc = 0;
74
mp_free(op->digs);
75
}
76
77
void
78
mpi_set(mpi *rop, mpi *op)
79
{
80
if (rop != op) {
81
if (rop->alloc < op->size) {
82
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op->size);
83
rop->alloc = op->size;
84
}
85
rop->size = op->size;
86
memcpy(rop->digs, op->digs, sizeof(BNS) * op->size);
87
rop->sign = op->sign;
88
}
89
}
90
91
void
92
mpi_seti(mpi *rop, long si)
93
{
94
unsigned long ui;
95
int sign = si < 0;
96
int size;
97
98
if (si == MINSLONG) {
99
ui = MINSLONG;
100
size = 2;
101
}
102
else {
103
if (sign)
104
ui = -si;
105
else
106
ui = si;
107
if (ui < CARRY)
108
size = 1;
109
else
110
size = 2;
111
}
112
113
if (rop->alloc < size) {
114
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size);
115
rop->alloc = size;
116
}
117
rop->size = size;
118
119
/* store data in small mp integer */
120
rop->digs[0] = (BNS)ui;
121
if (size > 1)
122
rop->digs[1] = (BNS)(ui >> BNSBITS);
123
rop->size = size;
124
125
/* adjust result sign */
126
rop->sign = sign;
127
}
128
129
static void
130
_mpi_seti(mpi *rop, long si)
131
{
132
unsigned long ui;
133
int sign = si < 0;
134
int size;
135
136
if (si == MINSLONG) {
137
ui = MINSLONG;
138
size = 2;
139
}
140
else {
141
if (sign)
142
ui = -si;
143
else
144
ui = si;
145
if (ui < CARRY)
146
size = 1;
147
else
148
size = 2;
149
}
150
151
rop->digs[0] = (BNS)ui;
152
if (size > 1)
153
rop->digs[1] = (BNS)(ui >> BNSBITS);
154
rop->size = size;
155
156
rop->sign = sign;
157
}
158
159
void
160
mpi_setd(mpi *rop, double d)
161
{
162
long i;
163
double mantissa;
164
int shift, exponent;
165
BNI size;
166
167
if (isnan(d))
168
d = 0.0;
169
else if (!finite(d))
170
d = copysign(1.0, d) * DBL_MAX;
171
172
/* check if number is larger than 1 */
173
if (fabs(d) < 1.0) {
174
rop->digs[0] = 0;
175
rop->size = 1;
176
rop->sign = d < 0.0;
177
178
return;
179
}
180
181
mantissa = frexp(d, &exponent);
182
if (mantissa < 0)
183
mantissa = -mantissa;
184
185
size = (exponent + (BNSBITS - 1)) / BNSBITS;
186
shift = BNSBITS - (exponent & (BNSBITS - 1));
187
188
/* adjust amount of memory */
189
if (rop->alloc < size) {
190
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size);
191
rop->alloc = size;
192
}
193
rop->size = size;
194
195
/* adjust the exponent */
196
if (shift < BNSBITS)
197
mantissa = ldexp(mantissa, -shift);
198
199
/* convert double */
200
for (i = size - 1; i >= 0 && mantissa != 0.0; i--) {
201
mantissa = ldexp(mantissa, BNSBITS);
202
rop->digs[i] = (BNS)mantissa;
203
mantissa -= rop->digs[i];
204
}
205
for (; i >= 0; i--)
206
rop->digs[i] = 0;
207
208
/* normalize */
209
if (size > 1 && rop->digs[size - 1] == 0)
210
--rop->size;
211
212
rop->sign = d < 0.0;
213
}
214
215
/* how many BNS in the given base, log(base) / log(CARRY) */
216
#ifdef LONG64
217
static double str_bases[37] = {
218
0.0000000000000000, 0.0000000000000000, 0.0312500000000000,
219
0.0495300781475362, 0.0625000000000000, 0.0725602529652301,
220
0.0807800781475362, 0.0877298413143002, 0.0937500000000000,
221
0.0990601562950723, 0.1038102529652301, 0.1081072380824156,
222
0.1120300781475362, 0.1156387411919092, 0.1189798413143002,
223
0.1220903311127662, 0.1250000000000000, 0.1277332137890731,
224
0.1303101562950723, 0.1327477347951120, 0.1350602529652300,
225
0.1372599194618363, 0.1393572380824156, 0.1413613111267817,
226
0.1432800781475362, 0.1451205059304602, 0.1468887411919092,
227
0.1485902344426084, 0.1502298413143002, 0.1518119060977367,
228
0.1533403311127662, 0.1548186346995899, 0.1562500000000000,
229
0.1576373162299517, 0.1589832137890731, 0.1602900942795302,
230
0.1615601562950723,
231
};
232
#else
233
static double str_bases[37] = {
234
0.0000000000000000, 0.0000000000000000, 0.0625000000000000,
235
0.0990601562950723, 0.1250000000000000, 0.1451205059304602,
236
0.1615601562950723, 0.1754596826286003, 0.1875000000000000,
237
0.1981203125901446, 0.2076205059304602, 0.2162144761648311,
238
0.2240601562950723, 0.2312774823838183, 0.2379596826286003,
239
0.2441806622255325, 0.2500000000000000, 0.2554664275781462,
240
0.2606203125901445, 0.2654954695902241, 0.2701205059304602,
241
0.2745198389236725, 0.2787144761648311, 0.2827226222535633,
242
0.2865601562950723, 0.2902410118609203, 0.2937774823838183,
243
0.2971804688852168, 0.3004596826286003, 0.3036238121954733,
244
0.3066806622255324, 0.3096372693991797, 0.3125000000000000,
245
0.3152746324599034, 0.3179664275781462, 0.3205801885590604,
246
0.3231203125901446,
247
};
248
#endif
249
250
void
251
mpi_setstr(mpi *rop, char *str, int base)
252
{
253
long i; /* counter */
254
int sign; /* result sign */
255
BNI carry; /* carry value */
256
BNI value; /* temporary value */
257
BNI size; /* size of result */
258
char *ptr; /* end of valid input */
259
260
/* initialization */
261
sign = 0;
262
carry = 0;
263
264
/* skip leading spaces */
265
while (isspace(*str))
266
++str;
267
268
/* check if sign supplied */
269
if (*str == '-') {
270
sign = 1;
271
++str;
272
}
273
else if (*str == '+')
274
++str;
275
276
/* skip leading zeros */
277
while (*str == '0')
278
++str;
279
280
ptr = str;
281
while (*ptr) {
282
if (*ptr >= '0' && *ptr <= '9') {
283
if (*ptr - '0' >= base)
284
break;
285
}
286
else if (*ptr >= 'A' && *ptr <= 'Z') {
287
if (*ptr - 'A' + 10 >= base)
288
break;
289
}
290
else if (*ptr >= 'a' && *ptr <= 'z') {
291
if (*ptr - 'a' + 10 >= base)
292
break;
293
}
294
else
295
break;
296
++ptr;
297
}
298
299
/* resulting size */
300
size = (ptr - str) * str_bases[base] + 1;
301
302
/* make sure rop has enough storage */
303
if (rop->alloc < size) {
304
rop->digs = mp_realloc(rop->digs, size * sizeof(BNS));
305
rop->alloc = size;
306
}
307
rop->size = size;
308
309
/* initialize rop to zero */
310
memset(rop->digs, '\0', size * sizeof(BNS));
311
312
/* set result sign */
313
rop->sign = sign;
314
315
/* convert string */
316
for (; str < ptr; str++) {
317
value = *str;
318
if (islower(value))
319
value = toupper(value);
320
value = value > '9' ? value - 'A' + 10 : value - '0';
321
value += rop->digs[0] * base;
322
carry = value >> BNSBITS;
323
rop->digs[0] = (BNS)value;
324
for (i = 1; i < size; i++) {
325
value = (BNI)rop->digs[i] * base + carry;
326
carry = value >> BNSBITS;
327
rop->digs[i] = (BNS)value;
328
}
329
}
330
331
/* normalize */
332
if (rop->size > 1 && rop->digs[rop->size - 1] == 0)
333
--rop->size;
334
}
335
336
void
337
mpi_add(mpi *rop, mpi *op1, mpi *op2)
338
{
339
mpi_addsub(rop, op1, op2, 0);
340
}
341
342
void
343
mpi_addi(mpi *rop, mpi *op1, long op2)
344
{
345
BNS digs[2];
346
mpi op;
347
348
op.digs = (BNS*)digs;
349
_mpi_seti(&op, op2);
350
351
mpi_addsub(rop, op1, &op, 0);
352
}
353
354
void
355
mpi_sub(mpi *rop, mpi *op1, mpi *op2)
356
{
357
mpi_addsub(rop, op1, op2, 1);
358
}
359
360
void
361
mpi_subi(mpi *rop, mpi *op1, long op2)
362
{
363
BNS digs[2];
364
mpi op;
365
366
op.digs = (BNS*)digs;
367
_mpi_seti(&op, op2);
368
369
mpi_addsub(rop, op1, &op, 1);
370
}
371
372
static void
373
mpi_addsub(mpi *rop, mpi *op1, mpi *op2, int sub)
374
{
375
long xlen; /* maximum result size */
376
377
if (sub ^ (op1->sign == op2->sign)) {
378
/* plus one for possible carry */
379
xlen = MAX(op1->size, op2->size) + 1;
380
if (rop->alloc < xlen) {
381
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xlen);
382
rop->alloc = xlen;
383
}
384
rop->size = mp_add(rop->digs, op1->digs, op2->digs,
385
op1->size, op2->size);
386
rop->sign = op1->sign;
387
}
388
else {
389
long cmp; /* check for larger operator */
390
391
cmp = mpi_cmpabs(op1, op2);
392
if (cmp == 0) {
393
rop->digs[0] = 0;
394
rop->size = 1;
395
rop->sign = 0;
396
}
397
else {
398
xlen = MAX(op1->size, op2->size);
399
if (rop->alloc < xlen) {
400
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xlen);
401
rop->alloc = xlen;
402
}
403
if (cmp > 0) {
404
rop->size = mp_sub(rop->digs, op1->digs, op2->digs,
405
op1->size, op2->size);
406
rop->sign = op1->sign;
407
}
408
else {
409
rop->size = mp_sub(rop->digs, op2->digs, op1->digs,
410
op2->size, op1->size);
411
rop->sign = sub ^ op2->sign;
412
}
413
}
414
}
415
}
416
417
void
418
mpi_mul(mpi *rop, mpi *op1, mpi *op2)
419
{
420
int sign; /* sign flag */
421
BNS *digs; /* result data */
422
long xsize; /* result size */
423
424
/* get result sign */
425
sign = op1->sign ^ op2->sign;
426
427
/* check for special cases */
428
if (op1->size == 1) {
429
if (*op1->digs == 0) {
430
/* multiply by 0 */
431
mpi_seti(rop, 0);
432
return;
433
}
434
else if (*op1->digs == 1) {
435
/* multiply by +-1 */
436
if (rop->alloc < op2->size) {
437
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op2->size);
438
rop->alloc = op2->size;
439
}
440
rop->size = op2->size;
441
memmove(rop->digs, op2->digs, sizeof(BNS) * op2->size);
442
rop->sign = op2->size > 1 || *op2->digs ? sign : 0;
443
444
return;
445
}
446
}
447
else if (op2->size == 1) {
448
if (*op2->digs == 0) {
449
/* multiply by 0 */
450
mpi_seti(rop, 0);
451
return;
452
}
453
else if (*op2->digs == 1) {
454
/* multiply by +-1 */
455
if (rop->alloc < op1->size) {
456
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op1->size);
457
rop->alloc = op1->size;
458
}
459
rop->size = op1->size;
460
memmove(rop->digs, op1->digs, sizeof(BNS) * op1->size);
461
rop->sign = op1->size > 1 || *op1->digs ? sign : 0;
462
463
return;
464
}
465
}
466
467
/* allocate result data and set it to zero */
468
xsize = op1->size + op2->size;
469
if (rop->digs == op1->digs || rop->digs == op2->digs)
470
/* rop is also an operand */
471
digs = mp_calloc(1, sizeof(BNS) * xsize);
472
else {
473
if (rop->alloc < xsize) {
474
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xsize);
475
rop->alloc = xsize;
476
}
477
digs = rop->digs;
478
memset(digs, '\0', sizeof(BNS) * xsize);
479
}
480
481
/* multiply operands */
482
xsize = mp_mul(digs, op1->digs, op2->digs, op1->size, op2->size);
483
484
/* store result in rop */
485
if (digs != rop->digs) {
486
/* if rop was an operand, free old data */
487
mp_free(rop->digs);
488
rop->digs = digs;
489
}
490
rop->size = xsize;
491
492
/* set result sign */
493
rop->sign = sign;
494
}
495
496
void
497
mpi_muli(mpi *rop, mpi *op1, long op2)
498
{
499
BNS digs[2];
500
mpi op;
501
502
op.digs = (BNS*)digs;
503
_mpi_seti(&op, op2);
504
505
mpi_mul(rop, op1, &op);
506
}
507
508
void
509
mpi_div(mpi *rop, mpi *num, mpi *den)
510
{
511
mpi_divqr(rop, NULL, num, den);
512
}
513
514
void
515
mpi_rem(mpi *rop, mpi *num, mpi *den)
516
{
517
mpi_divqr(NULL, rop, num, den);
518
}
519
520
/*
521
* Could/should be changed to not allocate qdigs if qrop is NULL
522
* Performance wouldn't suffer too much with a test on every loop iteration.
523
*/
524
void
525
mpi_divqr(mpi *qrop, mpi *rrop, mpi *num, mpi *den)
526
{
527
long i, j; /* counters */
528
int qsign; /* sign of quotient */
529
int rsign; /* sign of remainder */
530
BNI qsize; /* size of quotient */
531
BNI rsize; /* size of remainder */
532
BNS qest; /* estimative of quotient value */
533
BNS *qdigs, *rdigs; /* work copy or result */
534
BNS *ndigs, *ddigs; /* work copy or divisor and dividend */
535
BNI value; /* temporary result */
536
long svalue; /* signed temporary result (2's complement) */
537
BNS carry, scarry, denorm; /* carry and normalization */
538
BNI dpos, npos; /* offsets in data */
539
540
/* get signs */
541
rsign = num->sign;
542
qsign = rsign ^ den->sign;
543
544
/* check for special case */
545
if (num->size < den->size) {
546
/* quotient is zero and remainder is numerator */
547
if (rrop && rrop->digs != num->digs) {
548
if (rrop->alloc < num->size) {
549
rrop->digs = mp_realloc(rrop->digs, sizeof(BNS) * num->size);
550
rrop->alloc = num->size;
551
}
552
rrop->size = num->size;
553
memcpy(rrop->digs, num->digs, sizeof(BNS) * num->size);
554
rrop->sign = rsign;
555
}
556
if (qrop)
557
mpi_seti(qrop, 0);
558
559
return;
560
}
561
562
/* estimate result sizes */
563
rsize = den->size;
564
qsize = num->size - den->size + 1;
565
566
/* offsets */
567
npos = num->size - 1;
568
dpos = den->size - 1;
569
570
/* allocate space for quotient and remainder */
571
if (qrop == NULL || qrop->digs == num->digs || qrop->digs == den->digs)
572
qdigs = mp_calloc(1, sizeof(BNS) * qsize);
573
else {
574
if (qrop->alloc < qsize) {
575
qrop->digs = mp_realloc(qrop->digs, sizeof(BNS) * qsize);
576
qrop->alloc = qsize;
577
}
578
memset(qrop->digs, '\0', sizeof(BNS) * qsize);
579
qdigs = qrop->digs;
580
}
581
if (rrop) {
582
if (rrop->digs == num->digs || rrop->digs == den->digs)
583
rdigs = mp_calloc(1, sizeof(BNS) * rsize);
584
else {
585
if (rrop->alloc < rsize) {
586
rrop->digs = mp_realloc(rrop->digs, sizeof(BNS) * rsize);
587
rrop->alloc = rsize;
588
}
589
memset(rrop->digs, '\0', sizeof(BNS) * rsize);
590
rdigs = rrop->digs;
591
}
592
}
593
else
594
rdigs = NULL; /* fix gcc warning */
595
596
/* special case, only one word in divisor */
597
if (dpos == 0) {
598
for (carry = 0, i = npos; i >= 0; i--) {
599
value = ((BNI)carry << BNSBITS) + num->digs[i];
600
qdigs[i] = (BNS)(value / den->digs[0]);
601
carry = (BNS)(value % den->digs[0]);
602
}
603
if (rrop)
604
rdigs[0] = carry;
605
606
goto mpi_divqr_done;
607
}
608
609
/* make work copy of numerator */
610
ndigs = mp_malloc(sizeof(BNS) * (num->size + 1));
611
/* allocate one extra word an update offset */
612
memcpy(ndigs, num->digs, sizeof(BNS) * num->size);
613
ndigs[num->size] = 0;
614
++npos;
615
616
/* normalize */
617
denorm = (BNS)((BNI)CARRY / ((BNI)(den->digs[dpos]) + 1));
618
619
if (denorm > 1) {
620
/* i <= num->size because ndigs has an extra word */
621
for (carry = 0, i = 0; i <= num->size; i++) {
622
value = ndigs[i] * (BNI)denorm + carry;
623
ndigs[i] = (BNS)value;
624
carry = (BNS)(value >> BNSBITS);
625
}
626
/* make work copy of denominator */
627
ddigs = mp_malloc(sizeof(BNS) * den->size);
628
memcpy(ddigs, den->digs, sizeof(BNS) * den->size);
629
for (carry = 0, i = 0; i < den->size; i++) {
630
value = ddigs[i] * (BNI)denorm + carry;
631
ddigs[i] = (BNS)value;
632
carry = (BNS)(value >> BNSBITS);
633
}
634
}
635
else
636
/* only allocate copy of denominator if going to change it */
637
ddigs = den->digs;
638
639
/* divide mp integers */
640
for (j = qsize - 1; j >= 0; j--, npos--) {
641
/* estimate quotient */
642
if (ndigs[npos] == ddigs[dpos])
643
qest = (BNS)SMASK;
644
else
645
qest = (BNS)((((BNI)(ndigs[npos]) << 16) + ndigs[npos - 1]) /
646
ddigs[dpos]);
647
648
while ((value = ((BNI)(ndigs[npos]) << 16) + ndigs[npos - 1] -
649
qest * (BNI)(ddigs[dpos])) < CARRY &&
650
ddigs[dpos - 1] * (BNI)qest >
651
(value << BNSBITS) + ndigs[npos - 2])
652
--qest;
653
654
/* multiply and subtract */
655
carry = scarry = 0;
656
for (i = 0; i < den->size; i++) {
657
value = qest * (BNI)ddigs[i] + carry;
658
carry = (BNS)(value >> BNSBITS);
659
svalue = (long)ndigs[npos - dpos + i - 1] - (long)(value & SMASK) -
660
(long)scarry;
661
ndigs[npos - dpos + i - 1] = (BNS)svalue;
662
scarry = svalue < 0;
663
}
664
665
svalue = (long)ndigs[npos] - (long)(carry & SMASK) - (long)scarry;
666
ndigs[npos] = (BNS)svalue;
667
668
if (svalue & LMASK) {
669
/* quotient too big */
670
--qest;
671
carry = 0;
672
for (i = 0; i < den->size; i++) {
673
value = ndigs[npos - dpos + i - 1] + (BNI)carry + (BNI)ddigs[i];
674
ndigs[npos - dpos + i - 1] = (BNS)value;
675
carry = (BNS)(value >> BNSBITS);
676
}
677
ndigs[npos] += carry;
678
}
679
680
qdigs[j] = qest;
681
}
682
683
/* calculate remainder */
684
if (rrop) {
685
for (carry = 0, j = dpos; j >= 0; j--) {
686
value = ((BNI)carry << BNSBITS) + ndigs[j];
687
rdigs[j] = (BNS)(value / denorm);
688
carry = (BNS)(value % denorm);
689
}
690
}
691
692
mp_free(ndigs);
693
if (ddigs != den->digs)
694
mp_free(ddigs);
695
696
mpi_divqr_done:
697
if (rrop) {
698
if (rrop->digs != rdigs)
699
mp_free(rrop->digs);
700
/* normalize remainder */
701
for (i = rsize - 1; i >= 0; i--)
702
if (rdigs[i] != 0)
703
break;
704
if (i != rsize - 1) {
705
if (i < 0) {
706
rsign = 0;
707
rsize = 1;
708
}
709
else
710
rsize = i + 1;
711
}
712
rrop->digs = rdigs;
713
rrop->sign = rsign;
714
rrop->size = rsize;
715
}
716
717
/* normalize quotient */
718
if (qrop) {
719
if (qrop->digs != qdigs)
720
mp_free(qrop->digs);
721
for (i = qsize - 1; i >= 0; i--)
722
if (qdigs[i] != 0)
723
break;
724
if (i != qsize - 1) {
725
if (i < 0) {
726
qsign = 0;
727
qsize = 1;
728
}
729
else
730
qsize = i + 1;
731
}
732
qrop->digs = qdigs;
733
qrop->sign = qsign;
734
qrop->size = qsize;
735
}
736
else
737
mp_free(qdigs);
738
}
739
740
long
741
mpi_divqri(mpi *qrop, mpi *num, long den)
742
{
743
BNS ddigs[2];
744
mpi dop, rrop;
745
long remainder;
746
747
dop.digs = (BNS*)ddigs;
748
_mpi_seti(&dop, den);
749
750
memset(&rrop, '\0', sizeof(mpi));
751
mpi_init(&rrop);
752
mpi_divqr(qrop, &rrop, num, &dop);
753
remainder = rrop.digs[0];
754
if (rrop.size > 1)
755
remainder |= (BNI)(rrop.digs[1]) << BNSBITS;
756
if (rrop.sign)
757
remainder = -remainder;
758
mpi_clear(&rrop);
759
760
return (remainder);
761
}
762
763
void
764
mpi_divi(mpi *rop, mpi *num, long den)
765
{
766
BNS ddigs[2];
767
mpi dop;
768
769
dop.digs = (BNS*)ddigs;
770
_mpi_seti(&dop, den);
771
772
mpi_divqr(rop, NULL, num, &dop);
773
}
774
775
long
776
mpi_remi(mpi *num, long den)
777
{
778
return (mpi_divqri(NULL, num, den));
779
}
780
781
void
782
mpi_mod(mpi *rop, mpi *num, mpi *den)
783
{
784
mpi_rem(rop, num, den);
785
if (num->sign ^ den->sign)
786
mpi_add(rop, rop, den);
787
}
788
789
long
790
mpi_modi(mpi *num, long den)
791
{
792
long remainder;
793
794
remainder = mpi_remi(num, den);
795
if (num->sign ^ (den < 0))
796
remainder += den;
797
798
return (remainder);
799
}
800
801
void
802
mpi_gcd(mpi *rop, mpi *num, mpi *den)
803
{
804
long cmp;
805
mpi rem;
806
807
/* check if result already given */
808
cmp = mpi_cmpabs(num, den);
809
810
/* check if num is equal to den or if num is zero */
811
if (cmp == 0 || (num->size == 1 && num->digs[0] == 0)) {
812
mpi_set(rop, den);
813
rop->sign = 0;
814
return;
815
}
816
/* check if den is not zero */
817
if (den->size == 1 && den->digs[0] == 0) {
818
mpi_set(rop, num);
819
rop->sign = 0;
820
return;
821
}
822
823
/* don't call mpi_init, relies on realloc(0, size) == malloc(size) */
824
memset(&rem, '\0', sizeof(mpi));
825
826
/* if num larger than den */
827
if (cmp > 0) {
828
mpi_rem(&rem, num, den);
829
if (rem.size == 1 && rem.digs[0] == 0) {
830
/* exact division */
831
mpi_set(rop, den);
832
rop->sign = 0;
833
mpi_clear(&rem);
834
return;
835
}
836
mpi_set(rop, den);
837
}
838
else {
839
mpi_rem(&rem, den, num);
840
if (rem.size == 1 && rem.digs[0] == 0) {
841
/* exact division */
842
mpi_set(rop, num);
843
rop->sign = 0;
844
mpi_clear(&rem);
845
return;
846
}
847
mpi_set(rop, num);
848
}
849
850
/* loop using positive values */
851
rop->sign = rem.sign = 0;
852
853
/* cannot optimize this inverting rem/rop assignment earlier
854
* because rop mais be an operand */
855
mpi_swap(rop, &rem);
856
857
/* Euclides algorithm */
858
for (;;) {
859
mpi_rem(&rem, &rem, rop);
860
if (rem.size == 1 && rem.digs[0] == 0)
861
break;
862
mpi_swap(rop, &rem);
863
}
864
mpi_clear(&rem);
865
}
866
867
void
868
mpi_lcm(mpi *rop, mpi *num, mpi *den)
869
{
870
mpi gcd;
871
872
/* check for zero operand */
873
if ((num->size == 1 && num->digs[0] == 0) ||
874
(den->size == 1 && den->digs[0] == 0)) {
875
rop->digs[0] = 0;
876
rop->sign = 0;
877
return;
878
}
879
880
/* don't call mpi_init, relies on realloc(0, size) == malloc(size) */
881
memset(&gcd, '\0', sizeof(mpi));
882
883
mpi_gcd(&gcd, num, den);
884
mpi_div(&gcd, den, &gcd);
885
mpi_mul(rop, &gcd, num);
886
rop->sign = 0;
887
888
mpi_clear(&gcd);
889
}
890
891
void
892
mpi_pow(mpi *rop, mpi *op, unsigned long exp)
893
{
894
mpi zop, top;
895
896
if (exp == 2) {
897
mpi_mul(rop, op, op);
898
return;
899
}
900
/* check for op**0 */
901
else if (exp == 0) {
902
rop->digs[0] = 1;
903
rop->size = 1;
904
rop->sign = 0;
905
return;
906
}
907
else if (exp == 1) {
908
mpi_set(rop, op);
909
return;
910
}
911
else if (op->size == 1) {
912
if (op->digs[0] == 0) {
913
mpi_seti(rop, 0);
914
return;
915
}
916
else if (op->digs[0] == 1) {
917
mpi_seti(rop, op->sign && (exp & 1) ? -1 : 1);
918
return;
919
}
920
}
921
922
memset(&zop, '\0', sizeof(mpi));
923
memset(&top, '\0', sizeof(mpi));
924
mpi_set(&zop, op);
925
mpi_set(&top, op);
926
for (--exp; exp; exp >>= 1) {
927
if (exp & 1)
928
mpi_mul(&zop, &top, &zop);
929
mpi_mul(&top, &top, &top);
930
}
931
932
mpi_clear(&top);
933
rop->sign = zop.sign;
934
mp_free(rop->digs);
935
rop->digs = zop.digs;
936
rop->size = zop.size;
937
}
938
939
/* Find integer root of given number using the iteration
940
* x{n+1} = ((K-1) * x{n} + N / x{n}^(K-1)) / K
941
*/
942
int
943
mpi_root(mpi *rop, mpi *op, unsigned long nth)
944
{
945
long bits, cmp;
946
int exact;
947
int sign;
948
mpi *r, t, temp, quot, old, rem;
949
950
sign = op->sign;
951
952
/* divide by zero op**1/0 */
953
if (nth == 0) {
954
int one = 1, zero = 0;
955
one = one / zero;
956
}
957
/* result is complex */
958
if (sign && !(nth & 1)) {
959
int one = 1, zero = 0;
960
one = one / zero;
961
}
962
963
/* special case op**1/1 = op */
964
if (nth == 1) {
965
mpi_set(rop, op);
966
return (1);
967
}
968
969
bits = mpi_getsize(op, 2) - 2;
970
971
if (bits < 0 || bits / nth == 0) {
972
/* integral root is surely less than 2 */
973
exact = op->size == 1 && (op->digs[0] == 1 || op->digs[0] == 0);
974
mpi_seti(rop, sign ? -1 : op->digs[0] == 0 ? 0 : 1);
975
976
return (exact == 1);
977
}
978
979
/* initialize */
980
if (rop != op)
981
r = rop;
982
else {
983
r = &t;
984
memset(r, '\0', sizeof(mpi));
985
}
986
memset(&temp, '\0', sizeof(mpi));
987
memset(", '\0', sizeof(mpi));
988
memset(&old, '\0', sizeof(mpi));
989
memset(&rem, '\0', sizeof(mpi));
990
991
if (sign)
992
r->sign = 0;
993
994
/* root aproximation */
995
mpi_ash(r, op, -(bits - (bits / nth)));
996
997
for (;;) {
998
mpi_pow(&temp, r, nth - 1);
999
mpi_divqr(", &rem, op, &temp);
1000
cmp = mpi_cmpabs(r, ");
1001
if (cmp == 0) {
1002
exact = mpi_cmpi(&rem, 0) == 0;
1003
break;
1004
}
1005
else if (cmp < 0) {
1006
if (mpi_cmpabs(r, &old) == 0) {
1007
exact = 0;
1008
break;
1009
}
1010
mpi_set(&old, r);
1011
}
1012
mpi_muli(&temp, r, nth - 1);
1013
mpi_add(", ", &temp);
1014
mpi_divi(r, ", nth);
1015
}
1016
1017
mpi_clear(&temp);
1018
mpi_clear(");
1019
mpi_clear(&old);
1020
mpi_clear(&rem);
1021
if (r != rop) {
1022
mpi_set(rop, r);
1023
mpi_clear(r);
1024
}
1025
rop->sign = sign;
1026
1027
return (exact);
1028
}
1029
1030
/*
1031
* Find square root using the iteration:
1032
* x{n+1} = (x{n}+N/x{n})/2
1033
*/
1034
int
1035
mpi_sqrt(mpi *rop, mpi *op)
1036
{
1037
long bits, cmp;
1038
int exact;
1039
mpi *r, t, quot, rem, old;
1040
1041
/* result is complex */
1042
if (op->sign) {
1043
int one = 1, zero = 0;
1044
one = one / zero;
1045
}
1046
1047
bits = mpi_getsize(op, 2) - 2;
1048
1049
if (bits < 2) {
1050
/* integral root is surely less than 2 */
1051
exact = op->size == 1 && (op->digs[0] == 1 || op->digs[0] == 0);
1052
mpi_seti(rop, op->digs[0] == 0 ? 0 : 1);
1053
1054
return (exact == 1);
1055
}
1056
1057
/* initialize */
1058
if (rop != op)
1059
r = rop;
1060
else {
1061
r = &t;
1062
memset(r, '\0', sizeof(mpi));
1063
}
1064
memset(", '\0', sizeof(mpi));
1065
memset(&rem, '\0', sizeof(mpi));
1066
memset(&old, '\0', sizeof(mpi));
1067
1068
/* root aproximation */
1069
mpi_ash(r, op, -(bits - (bits / 2)));
1070
1071
for (;;) {
1072
if (mpi_cmpabs(r, &old) == 0) {
1073
exact = 0;
1074
break;
1075
}
1076
mpi_divqr(", &rem, op, r);
1077
cmp = mpi_cmpabs(", r);
1078
if (cmp == 0) {
1079
exact = mpi_cmpi(&rem, 0) == 0;
1080
break;
1081
}
1082
else if (cmp > 0 && rem.size == 1 && rem.digs[0] == 0)
1083
mpi_subi(", ", 1);
1084
mpi_set(&old, r);
1085
mpi_add(r, r, ");
1086
mpi_ash(r, r, -1);
1087
}
1088
mpi_clear(");
1089
mpi_clear(&rem);
1090
mpi_clear(&old);
1091
if (r != rop) {
1092
mpi_set(rop, r);
1093
mpi_clear(r);
1094
}
1095
1096
return (exact);
1097
}
1098
1099
void
1100
mpi_ash(mpi *rop, mpi *op, long shift)
1101
{
1102
long i; /* counter */
1103
long xsize; /* maximum result size */
1104
BNS *digs;
1105
1106
/* check for 0 shift, multiply/divide by 1 */
1107
if (shift == 0) {
1108
if (rop != op) {
1109
if (rop->alloc < op->size) {
1110
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op->size);
1111
rop->alloc = op->size;
1112
}
1113
rop->size = op->size;
1114
memcpy(rop->digs, op->digs, sizeof(BNS) * op->size);
1115
}
1116
1117
return;
1118
}
1119
else if (op->size == 1 && op->digs[0] == 0) {
1120
rop->sign = 0;
1121
rop->size = 1;
1122
rop->digs[0] = 0;
1123
1124
return;
1125
}
1126
1127
/* check shift and initialize */
1128
if (shift > 0)
1129
xsize = op->size + (shift / BNSBITS) + 1;
1130
else {
1131
xsize = op->size - ((-shift) / BNSBITS);
1132
if (xsize <= 0) {
1133
rop->size = 1;
1134
rop->sign = op->sign;
1135
rop->digs[0] = op->sign ? 1 : 0;
1136
1137
return;
1138
}
1139
}
1140
1141
/* allocate/adjust memory for result */
1142
if (rop == op)
1143
digs = mp_malloc(sizeof(BNS) * xsize);
1144
else {
1145
if (rop->alloc < xsize) {
1146
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xsize);
1147
rop->alloc = xsize;
1148
}
1149
digs = rop->digs;
1150
}
1151
1152
/* left shift, multiply by power of two */
1153
if (shift > 0)
1154
rop->size = mp_lshift(digs, op->digs, op->size, shift);
1155
/* right shift, divide by power of two */
1156
else {
1157
long carry = 0;
1158
1159
if (op->sign) {
1160
BNI words, bits;
1161
1162
words = -shift / BNSBITS;
1163
bits = -shift % BNSBITS;
1164
for (i = 0; i < words; i++)
1165
carry |= op->digs[xsize + i];
1166
if (!carry) {
1167
for (i = 0; i < bits; i++)
1168
if (op->digs[op->size - xsize] & (1 << i)) {
1169
carry = 1;
1170
break;
1171
}
1172
}
1173
}
1174
rop->size = mp_rshift(digs, op->digs, op->size, -shift);
1175
1176
if (carry)
1177
/* emulates two's complement subtracting 1 from the result */
1178
rop->size = mp_add(digs, digs, mpone.digs, rop->size, 1);
1179
}
1180
1181
if (rop->digs != digs) {
1182
mp_free(rop->digs);
1183
rop->alloc = rop->size;
1184
rop->digs = digs;
1185
}
1186
rop->sign = op->sign;
1187
}
1188
1189
static INLINE BNS
1190
mpi_logic(BNS op1, BNS op2, BNS op)
1191
{
1192
switch (op) {
1193
case '&':
1194
return (op1 & op2);
1195
case '|':
1196
return (op1 | op2);
1197
case '^':
1198
return (op1 ^ op2);
1199
}
1200
1201
return (SMASK);
1202
}
1203
1204
static void
1205
mpi_log(mpi *rop, mpi *op1, mpi *op2, BNS op)
1206
{
1207
long i; /* counter */
1208
long c, c1, c2; /* carry */
1209
BNS *digs, *digs1, *digs2; /* pointers to mp data */
1210
BNI size, size1, size2;
1211
BNS sign, sign1, sign2;
1212
BNS n, n1, n2; /* logical operands */
1213
BNI sum;
1214
1215
/* initialize */
1216
size1 = op1->size;
1217
size2 = op2->size;
1218
1219
sign1 = op1->sign ? SMASK : 0;
1220
sign2 = op2->sign ? SMASK : 0;
1221
1222
sign = mpi_logic(sign1, sign2, op);
1223
1224
size = MAX(size1, size2);
1225
if (sign)
1226
++size;
1227
if (rop->alloc < size) {
1228
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size);
1229
rop->alloc = size;
1230
}
1231
1232
digs = rop->digs;
1233
digs1 = op1->digs;
1234
digs2 = op2->digs;
1235
1236
c = c1 = c2 = 1;
1237
1238
/* apply logical operation */
1239
for (i = 0; i < size; i++) {
1240
if (i >= size1)
1241
n1 = sign1;
1242
else if (sign1) {
1243
sum = (BNI)(BNS)(~digs1[i]) + c1;
1244
c1 = (long)(sum >> BNSBITS);
1245
n1 = (BNS)sum;
1246
}
1247
else
1248
n1 = digs1[i];
1249
1250
if (i >= size2)
1251
n2 = sign2;
1252
else if (sign2) {
1253
sum = (BNI)(BNS)(~digs2[i]) + c2;
1254
c2 = (long)(sum >> BNSBITS);
1255
n2 = (BNS)sum;
1256
}
1257
else
1258
n2 = digs2[i];
1259
1260
n = mpi_logic(n1, n2, op);
1261
if (sign) {
1262
sum = (BNI)(BNS)(~n) + c;
1263
c = (long)(sum >> BNSBITS);
1264
digs[i] = (BNS)sum;
1265
}
1266
else
1267
digs[i] = n;
1268
}
1269
1270
/* normalize */
1271
for (i = size - 1; i >= 0; i--)
1272
if (digs[i] != 0)
1273
break;
1274
if (i != size - 1) {
1275
if (i < 0) {
1276
sign = 0;
1277
size = 1;
1278
}
1279
else
1280
size = i + 1;
1281
}
1282
1283
rop->sign = sign != 0;
1284
rop->size = size;
1285
}
1286
1287
void
1288
mpi_and(mpi *rop, mpi *op1, mpi *op2)
1289
{
1290
mpi_log(rop, op1, op2, '&');
1291
}
1292
1293
void
1294
mpi_ior(mpi *rop, mpi *op1, mpi *op2)
1295
{
1296
mpi_log(rop, op1, op2, '|');
1297
}
1298
1299
void
1300
mpi_xor(mpi *rop, mpi *op1, mpi *op2)
1301
{
1302
mpi_log(rop, op1, op2, '^');
1303
}
1304
1305
void
1306
mpi_com(mpi *rop, mpi *op)
1307
{
1308
static BNS digs[1] = { 1 };
1309
static mpi one = { 0, 1, 1, (BNS*)&digs };
1310
1311
mpi_log(rop, rop, &one, '^');
1312
}
1313
1314
void
1315
mpi_neg(mpi *rop, mpi *op)
1316
{
1317
int sign = op->sign ^ 1;
1318
1319
if (rop->digs != op->digs) {
1320
if (rop->alloc < op->size) {
1321
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * rop->size);
1322
rop->alloc = op->size;
1323
}
1324
rop->size = op->size;
1325
memcpy(rop->digs, op->digs, sizeof(BNS) * rop->size);
1326
}
1327
1328
rop->sign = sign;
1329
}
1330
1331
void
1332
mpi_abs(mpi *rop, mpi *op)
1333
{
1334
if (rop->digs != op->digs) {
1335
if (rop->alloc < op->size) {
1336
rop->digs = mp_realloc(rop->digs, sizeof(BNS) * rop->size);
1337
rop->alloc = op->size;
1338
}
1339
rop->size = op->size;
1340
memcpy(rop->digs, op->digs, sizeof(BNS) * rop->size);
1341
}
1342
1343
rop->sign = 0;
1344
}
1345
1346
int
1347
mpi_cmp(mpi *op1, mpi *op2)
1348
{
1349
if (op1->sign ^ op2->sign)
1350
return (op1->sign ? -1 : 1);
1351
1352
if (op1->size == op2->size) {
1353
long i, cmp = 0;
1354
1355
for (i = op1->size - 1; i >= 0; i--)
1356
if ((cmp = (long)op1->digs[i] - (long)op2->digs[i]) != 0)
1357
break;
1358
1359
return (cmp == 0 ? 0 : (cmp < 0) ^ op1->sign ? -1 : 1);
1360
}
1361
1362
return ((op1->size < op2->size) ^ op1->sign ? -1 : 1);
1363
}
1364
1365
int
1366
mpi_cmpi(mpi *op1, long op2)
1367
{
1368
long cmp;
1369
1370
if (op1->size > 2)
1371
return (op1->sign ? -1 : 1);
1372
1373
cmp = op1->digs[0];
1374
if (op1->size == 2) {
1375
cmp |= (long)op1->digs[1] << BNSBITS;
1376
if (cmp == MINSLONG)
1377
return (op2 == cmp && op1->sign ? 0 : op1->sign ? -1 : 1);
1378
}
1379
if (op1->sign)
1380
cmp = -cmp;
1381
1382
return (cmp - op2);
1383
}
1384
1385
int
1386
mpi_cmpabs(mpi *op1, mpi *op2)
1387
{
1388
if (op1->size == op2->size) {
1389
long i, cmp = 0;
1390
1391
for (i = op1->size - 1; i >= 0; i--)
1392
if ((cmp = (long)op1->digs[i] - (long)op2->digs[i]) != 0)
1393
break;
1394
1395
return (cmp);
1396
}
1397
1398
return ((op1->size < op2->size) ? -1 : 1);
1399
}
1400
1401
int
1402
mpi_cmpabsi(mpi *op1, long op2)
1403
{
1404
unsigned long cmp;
1405
1406
if (op1->size > 2)
1407
return (1);
1408
1409
cmp = op1->digs[0];
1410
if (op1->size == 2)
1411
cmp |= (unsigned long)op1->digs[1] << BNSBITS;
1412
1413
return (cmp > op2 ? 1 : cmp == op2 ? 0 : -1);
1414
}
1415
1416
int
1417
mpi_sgn(mpi *op)
1418
{
1419
return (op->sign ? -1 : op->size > 1 || op->digs[0] ? 1 : 0);
1420
}
1421
1422
void
1423
mpi_swap(mpi *op1, mpi *op2)
1424
{
1425
if (op1 != op2) {
1426
mpi swap;
1427
1428
memcpy(&swap, op1, sizeof(mpi));
1429
memcpy(op1, op2, sizeof(mpi));
1430
memcpy(op2, &swap, sizeof(mpi));
1431
}
1432
}
1433
1434
int
1435
mpi_fiti(mpi *op)
1436
{
1437
if (op->size == 1)
1438
return (1);
1439
else if (op->size == 2) {
1440
unsigned long value = ((BNI)(op->digs[1]) << BNSBITS) | op->digs[0];
1441
1442
if (value & MINSLONG)
1443
return (op->sign && value == MINSLONG) ? 1 : 0;
1444
1445
return (1);
1446
}
1447
1448
return (0);
1449
}
1450
1451
long
1452
mpi_geti(mpi *op)
1453
{
1454
long value;
1455
1456
value = op->digs[0];
1457
if (op->size > 1)
1458
value |= (BNI)(op->digs[1]) << BNSBITS;
1459
1460
return (op->sign && value != MINSLONG ? -value : value);
1461
}
1462
1463
double
1464
mpi_getd(mpi *op)
1465
{
1466
long i, len;
1467
double d = 0.0;
1468
int exponent;
1469
1470
#define FLOATDIGS sizeof(double) / sizeof(BNS)
1471
1472
switch (op->size) {
1473
case 2:
1474
d = (BNI)(op->digs[1]) << BNSBITS;
1475
case 1:
1476
d += op->digs[0];
1477
return (op->sign ? -d : d);
1478
default:
1479
break;
1480
}
1481
1482
for (i = 0, len = op->size; len > 0 && i < FLOATDIGS; i++)
1483
d = ldexp(d, BNSBITS) + op->digs[--len];
1484
d = frexp(d, &exponent);
1485
if (len > 0)
1486
exponent += len * BNSBITS;
1487
1488
if (d == 0.0)
1489
return (d);
1490
1491
d = ldexp(d, exponent);
1492
1493
return (op->sign ? -d : d);
1494
}
1495
1496
/* how many digits in a given base, floor(log(CARRY) / log(base)) */
1497
#ifdef LONG64
1498
static char dig_bases[37] = {
1499
0, 0, 32, 20, 16, 13, 12, 11, 10, 10, 9, 9, 8, 8, 8, 8,
1500
8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6,
1501
6, 6, 6, 6, 6,
1502
};
1503
#else
1504
static char dig_bases[37] = {
1505
0, 0, 16, 10, 8, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4,
1506
4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
1507
3, 3, 3, 3, 3,
1508
};
1509
#endif
1510
1511
/* how many digits per bit in a given base, log(2) / log(base) */
1512
static double bit_bases[37] = {
1513
0.0000000000000000, 0.0000000000000000, 1.0000000000000000,
1514
0.6309297535714575, 0.5000000000000000, 0.4306765580733931,
1515
0.3868528072345416, 0.3562071871080222, 0.3333333333333334,
1516
0.3154648767857287, 0.3010299956639811, 0.2890648263178878,
1517
0.2789429456511298, 0.2702381544273197, 0.2626495350371936,
1518
0.2559580248098155, 0.2500000000000000, 0.2446505421182260,
1519
0.2398124665681315, 0.2354089133666382, 0.2313782131597592,
1520
0.2276702486969530, 0.2242438242175754, 0.2210647294575037,
1521
0.2181042919855316, 0.2153382790366965, 0.2127460535533632,
1522
0.2103099178571525, 0.2080145976765095, 0.2058468324604344,
1523
0.2037950470905062, 0.2018490865820999, 0.2000000000000000,
1524
0.1982398631705605, 0.1965616322328226, 0.1949590218937863,
1525
0.1934264036172708,
1526
};
1527
1528
/* normalization base for string conversion, pow(base, dig_bases[base]) & ~CARRY */
1529
#ifdef LONG64
1530
static BNS big_bases[37] = {
1531
0x00000001, 0x00000001, 0x00000000, 0xCFD41B91, 0x00000000, 0x48C27395,
1532
0x81BF1000, 0x75DB9C97, 0x40000000, 0xCFD41B91, 0x3B9ACA00, 0x8C8B6D2B,
1533
0x19A10000, 0x309F1021, 0x57F6C100, 0x98C29B81, 0x00000000, 0x18754571,
1534
0x247DBC80, 0x3547667B, 0x4C4B4000, 0x6B5A6E1D, 0x94ACE180, 0xCAF18367,
1535
0x0B640000, 0x0E8D4A51, 0x1269AE40, 0x17179149, 0x1CB91000, 0x23744899,
1536
0x2B73A840, 0x34E63B41, 0x40000000, 0x4CFA3CC1, 0x5C13D840, 0x6D91B519,
1537
0x81BF1000,
1538
};
1539
#else
1540
static BNS big_bases[37] = {
1541
0x0001, 0x0001, 0x0000, 0xE6A9, 0x0000, 0x3D09, 0xB640, 0x41A7, 0x8000,
1542
0xE6A9, 0x2710, 0x3931, 0x5100, 0x6F91, 0x9610, 0xC5C1, 0x0000, 0x1331,
1543
0x16C8, 0x1ACB, 0x1F40, 0x242D, 0x2998, 0x2F87, 0x3600, 0x3D09, 0x44A8,
1544
0x4CE3, 0x55C0, 0x5F45, 0x6978, 0x745F, 0x8000, 0x8C61, 0x9988, 0xA77B,
1545
0xb640,
1546
};
1547
#endif
1548
1549
unsigned long
1550
mpi_getsize(mpi *op, int base)
1551
{
1552
unsigned long value, bits;
1553
1554
value = op->digs[op->size - 1];
1555
1556
/* count leading bits */
1557
if (value) {
1558
unsigned long count, carry;
1559
1560
for (count = 0, carry = CARRY >> 1; carry; count++, carry >>= 1)
1561
if (value & carry)
1562
break;
1563
1564
bits = BNSBITS - count;
1565
}
1566
else
1567
bits = 0;
1568
1569
return ((bits + (op->size - 1) * BNSBITS) * bit_bases[base] + 1);
1570
}
1571
1572
char *
1573
mpi_getstr(char *str, mpi *op, int base)
1574
{
1575
long i; /* counter */
1576
BNS *digs, *xdigs; /* copy of op data */
1577
BNI size; /* size of op */
1578
BNI digits; /* digits per word in given base */
1579
BNI bigbase; /* big base of given base */
1580
BNI strsize; /* size of resulting string */
1581
char *cp; /* pointer in str for conversion */
1582
1583
/* initialize */
1584
size = op->size;
1585
strsize = mpi_getsize(op, base) + op->sign + 1;
1586
1587
if (str == NULL)
1588
str = mp_malloc(strsize);
1589
1590
/* check for zero */
1591
if (size == 1 && op->digs[0] == 0) {
1592
str[0] = '0';
1593
str[1] = '\0';
1594
1595
return (str);
1596
}
1597
1598
digits = dig_bases[base];
1599
bigbase = big_bases[base];
1600
1601
cp = str + strsize;
1602
*--cp = '\0';
1603
1604
/* make copy of op data and adjust digs */
1605
xdigs = mp_malloc(size * sizeof(BNS));
1606
memcpy(xdigs, op->digs, size * sizeof(unsigned short));
1607
digs = xdigs + size - 1;
1608
1609
/* convert to string */
1610
for (;;) {
1611
long count = -1;
1612
BNI value;
1613
BNS quotient, remainder = 0;
1614
1615
/* if power of two base */
1616
if ((base & (base - 1)) == 0) {
1617
for (i = 0; i < size; i++) {
1618
quotient = remainder;
1619
remainder = digs[-i];
1620
digs[-i] = quotient;
1621
if (count < 0 && quotient)
1622
count = i;
1623
}
1624
}
1625
else {
1626
for (i = 0; i < size; i++) {
1627
value = digs[-i] + ((BNI)remainder << BNSBITS);
1628
quotient = (BNS)(value / bigbase);
1629
remainder = (BNS)(value % bigbase);
1630
digs[-i] = quotient;
1631
if (count < 0 && quotient)
1632
count = i;
1633
}
1634
}
1635
quotient = remainder;
1636
for (i = 0; i < digits; i++) {
1637
if (quotient == 0 && count < 0)
1638
break;
1639
remainder = quotient % base;
1640
quotient /= base;
1641
*--cp = remainder < 10 ? remainder + '0' : remainder - 10 + 'A';
1642
}
1643
if (count < 0)
1644
break;
1645
digs -= count;
1646
size -= count;
1647
}
1648
1649
/* adjust sign */
1650
if (op->sign)
1651
*--cp = '-';
1652
1653
/* remove any extra characters */
1654
if (cp > str)
1655
strcpy(str, cp);
1656
1657
mp_free(xdigs);
1658
1659
return (str);
1660
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1