1 // Copyright 2009 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
4
5 #include "go.h"
6
7 //
8 // return the significant
9 // words of the argument
10 //
11 static int
12 mplen(Mpint *a)
13 {
14 int i, n;
15 long *a1;
16
17 n = -1;
18 a1 = &a->a[0];
19 for(i=0; i<Mpprec; i++) {
20 if(*a1++ != 0)
21 n = i;
22 }
23 return n+1;
24 }
25
26 //
27 // left shift mpint by one
28 // ignores sign and overflow
29 //
30 static void
31 mplsh(Mpint *a)
32 {
33 long *a1, x;
34 int i, c;
35
36 c = 0;
37 a1 = &a->a[0];
38 for(i=0; i<Mpprec; i++) {
39 x = (*a1 << 1) + c;
40 c = 0;
41 if(x >= Mpbase) {
42 x -= Mpbase;
43 c = 1;
44 }
45 *a1++ = x;
46 }
47 }
48
49 //
50 // left shift mpint by Mpscale
51 // ignores sign and overflow
52 //
53 static void
54 mplshw(Mpint *a)
55 {
56 long *a1;
57 int i;
58
59 a1 = &a->a[Mpprec-1];
60 for(i=1; i<Mpprec; i++) {
61 a1[0] = a1[-1];
62 a1--;
63 }
64 a1[0] = 0;
65 }
66
67 //
68 // right shift mpint by one
69 // ignores sign and overflow
70 //
71 static void
72 mprsh(Mpint *a)
73 {
74 long *a1, x, lo;
75 int i, c;
76
77 c = 0;
78 lo = a->a[0] & 1;
79 a1 = &a->a[Mpprec];
80 for(i=0; i<Mpprec; i++) {
81 x = *--a1;
82 *a1 = (x + c) >> 1;
83 c = 0;
84 if(x & 1)
85 c = Mpbase;
86 }
87 if(a->neg && lo != 0)
88 mpaddcfix(a, -1);
89 }
90
91 //
92 // right shift mpint by Mpscale
93 // ignores sign and overflow
94 //
95 static void
96 mprshw(Mpint *a)
97 {
98 long *a1, lo;
99 int i;
100
101 lo = a->a[0];
102 a1 = &a->a[0];
103 for(i=1; i<Mpprec; i++) {
104 a1[0] = a1[1];
105 a1++;
106 }
107 a1[0] = 0;
108 if(a->neg && lo != 0)
109 mpaddcfix(a, -1);
110 }
111
112 //
113 // return the sign of (abs(a)-abs(b))
114 //
115 static int
116 mpcmp(Mpint *a, Mpint *b)
117 {
118 long x, *a1, *b1;
119 int i;
120
121 if(a->ovf || b->ovf) {
122 yyerror("ovf in cmp");
123 return 0;
124 }
125
126 a1 = &a->a[0] + Mpprec;
127 b1 = &b->a[0] + Mpprec;
128
129 for(i=0; i<Mpprec; i++) {
130 x = *--a1 - *--b1;
131 if(x > 0)
132 return +1;
133 if(x < 0)
134 return -1;
135 }
136 return 0;
137 }
138
139 //
140 // negate a
141 // ignore sign and ovf
142 //
143 static void
144 mpneg(Mpint *a)
145 {
146 long x, *a1;
147 int i, c;
148
149 a1 = &a->a[0];
150 c = 0;
151 for(i=0; i<Mpprec; i++) {
152 x = -*a1 -c;
153 c = 0;
154 if(x < 0) {
155 x += Mpbase;
156 c = 1;
157 }
158 *a1++ = x;
159 }
160 }
161
162 // shift left by s (or right by -s)
163 void
164 mpshiftfix(Mpint *a, int s)
165 {
166 if(s >= 0) {
167 while(s >= Mpscale) {
168 mplshw(a);
169 s -= Mpscale;
170 }
171 while(s > 0) {
172 mplsh(a);
173 s--;
174 }
175 } else {
176 s = -s;
177 while(s >= Mpscale) {
178 mprshw(a);
179 s -= Mpscale;
180 }
181 while(s > 0) {
182 mprsh(a);
183 s--;
184 }
185 }
186 }
187
188 /// implements fix arihmetic
189
190 void
191 mpaddfixfix(Mpint *a, Mpint *b)
192 {
193 int i, c;
194 long x, *a1, *b1;
195
196 if(a->ovf || b->ovf) {
197 yyerror("ovf in mpaddxx");
198 a->ovf = 1;
199 return;
200 }
201
202 c = 0;
203 a1 = &a->a[0];
204 b1 = &b->a[0];
205 if(a->neg != b->neg)
206 goto sub;
207
208 // perform a+b
209 for(i=0; i<Mpprec; i++) {
210 x = *a1 + *b1++ + c;
211 c = 0;
212 if(x >= Mpbase) {
213 x -= Mpbase;
214 c = 1;
215 }
216 *a1++ = x;
217 }
218 a->ovf = c;
219 if(a->ovf)
220 yyerror("set ovf in mpaddxx");
221
222 return;
223
224 sub:
225 // perform a-b
226 switch(mpcmp(a, b)) {
227 case 0:
228 mpmovecfix(a, 0);
229 break;
230
231 case 1:
232 for(i=0; i<Mpprec; i++) {
233 x = *a1 - *b1++ - c;
234 c = 0;
235 if(x < 0) {
236 x += Mpbase;
237 c = 1;
238 }
239 *a1++ = x;
240 }
241 break;
242
243 case -1:
244 a->neg ^= 1;
245 for(i=0; i<Mpprec; i++) {
246 x = *b1++ - *a1 - c;
247 c = 0;
248 if(x < 0) {
249 x += Mpbase;
250 c = 1;
251 }
252 *a1++ = x;
253 }
254 break;
255 }
256 }
257
258 void
259 mpmulfixfix(Mpint *a, Mpint *b)
260 {
261
262 int i, j, na, nb;
263 long *a1, x;
264 Mpint s, q;
265
266 if(a->ovf || b->ovf) {
267 yyerror("ovf in mpmulfixfix");
268 a->ovf = 1;
269 return;
270 }
271
272 // pick the smaller
273 // to test for bits
274 na = mplen(a);
275 nb = mplen(b);
276 if(na > nb) {
277 mpmovefixfix(&s, a);
278 a1 = &b->a[0];
279 na = nb;
280 } else {
281 mpmovefixfix(&s, b);
282 a1 = &a->a[0];
283 }
284 s.neg = 0;
285
286 mpmovecfix(&q, 0);
287 for(i=0; i<na; i++) {
288 x = *a1++;
289 for(j=0; j<Mpscale; j++) {
290 if(x & 1)
291 mpaddfixfix(&q, &s);
292 mplsh(&s);
293 x >>= 1;
294 }
295 }
296
297 q.neg = a->neg ^ b->neg;
298 mpmovefixfix(a, &q);
299 if(a->ovf)
300 yyerror("set ovf in mpmulfixfix");
301 }
302
303 void
304 mpmulfract(Mpint *a, Mpint *b)
305 {
306
307 int i, j;
308 long *a1, x;
309 Mpint s, q;
310
311 if(a->ovf || b->ovf) {
312 yyerror("ovf in mpmulflt");
313 a->ovf = 1;
314 return;
315 }
316
317 mpmovefixfix(&s, b);
318 a1 = &a->a[Mpprec];
319 s.neg = 0;
320 mpmovecfix(&q, 0);
321
322 x = *--a1;
323 if(x != 0)
324 yyerror("mpmulfract not normal");
325
326 for(i=0; i<Mpprec-1; i++) {
327 x = *--a1;
328 if(x == 0) {
329 mprshw(&s);
330 continue;
331 }
332 for(j=0; j<Mpscale; j++) {
333 x <<= 1;
334 if(x & Mpbase)
335 mpaddfixfix(&q, &s);
336 mprsh(&s);
337 }
338 }
339
340 q.neg = a->neg ^ b->neg;
341 mpmovefixfix(a, &q);
342 if(a->ovf)
343 yyerror("set ovf in mpmulflt");
344 }
345
346 void
347 mporfixfix(Mpint *a, Mpint *b)
348 {
349 int i;
350 long x, *a1, *b1;
351
352 if(a->ovf || b->ovf) {
353 yyerror("ovf in mporfixfix");
354 mpmovecfix(a, 0);
355 a->ovf = 1;
356 return;
357 }
358 if(a->neg) {
359 a->neg = 0;
360 mpneg(a);
361 }
362 if(b->neg)
363 mpneg(b);
364
365 a1 = &a->a[0];
366 b1 = &b->a[0];
367 for(i=0; i<Mpprec; i++) {
368 x = *a1 | *b1++;
369 *a1++ = x;
370 }
371
372 if(b->neg)
373 mpneg(b);
374 if(x & Mpsign) {
375 a->neg = 1;
376 mpneg(a);
377 }
378 }
379
380 void
381 mpandfixfix(Mpint *a, Mpint *b)
382 {
383 int i;
384 long x, *a1, *b1;
385
386 if(a->ovf || b->ovf) {
387 yyerror("ovf in mpandfixfix");
388 mpmovecfix(a, 0);
389 a->ovf = 1;
390 return;
391 }
392 if(a->neg) {
393 a->neg = 0;
394 mpneg(a);
395 }
396 if(b->neg)
397 mpneg(b);
398
399 a1 = &a->a[0];
400 b1 = &b->a[0];
401 for(i=0; i<Mpprec; i++) {
402 x = *a1 & *b1++;
403 *a1++ = x;
404 }
405
406 if(b->neg)
407 mpneg(b);
408 if(x & Mpsign) {
409 a->neg = 1;
410 mpneg(a);
411 }
412 }
413
414 void
415 mpandnotfixfix(Mpint *a, Mpint *b)
416 {
417 int i;
418 long x, *a1, *b1;
419
420 if(a->ovf || b->ovf) {
421 yyerror("ovf in mpandnotfixfix");
422 mpmovecfix(a, 0);
423 a->ovf = 1;
424 return;
425 }
426 if(a->neg) {
427 a->neg = 0;
428 mpneg(a);
429 }
430 if(b->neg)
431 mpneg(b);
432
433 a1 = &a->a[0];
434 b1 = &b->a[0];
435 for(i=0; i<Mpprec; i++) {
436 x = *a1 & ~*b1++;
437 *a1++ = x;
438 }
439
440 if(b->neg)
441 mpneg(b);
442 if(x & Mpsign) {
443 a->neg = 1;
444 mpneg(a);
445 }
446 }
447
448 void
449 mpxorfixfix(Mpint *a, Mpint *b)
450 {
451 int i;
452 long x, *a1, *b1;
453
454 if(a->ovf || b->ovf) {
455 yyerror("ovf in mporfixfix");
456 mpmovecfix(a, 0);
457 a->ovf = 1;
458 return;
459 }
460 if(a->neg) {
461 a->neg = 0;
462 mpneg(a);
463 }
464 if(b->neg)
465 mpneg(b);
466
467 a1 = &a->a[0];
468 b1 = &b->a[0];
469 for(i=0; i<Mpprec; i++) {
470 x = *a1 ^ *b1++;
471 *a1++ = x;
472 }
473
474 if(b->neg)
475 mpneg(b);
476 if(x & Mpsign) {
477 a->neg = 1;
478 mpneg(a);
479 }
480 }
481
482 void
483 mplshfixfix(Mpint *a, Mpint *b)
484 {
485 vlong s;
486
487 if(a->ovf || b->ovf) {
488 yyerror("ovf in mporfixfix");
489 mpmovecfix(a, 0);
490 a->ovf = 1;
491 return;
492 }
493 s = mpgetfix(b);
494 if(s < 0 || s >= Mpprec*Mpscale) {
495 yyerror("stupid shift: %lld", s);
496 mpmovecfix(a, 0);
497 return;
498 }
499
500 mpshiftfix(a, s);
501 }
502
503 void
504 mprshfixfix(Mpint *a, Mpint *b)
505 {
506 vlong s;
507
508 if(a->ovf || b->ovf) {
509 yyerror("ovf in mprshfixfix");
510 mpmovecfix(a, 0);
511 a->ovf = 1;
512 return;
513 }
514 s = mpgetfix(b);
515 if(s < 0 || s >= Mpprec*Mpscale) {
516 yyerror("stupid shift: %lld", s);
517 if(a->neg)
518 mpmovecfix(a, -1);
519 else
520 mpmovecfix(a, 0);
521 return;
522 }
523
524 mpshiftfix(a, -s);
525 }
526
527 void
528 mpnegfix(Mpint *a)
529 {
530 a->neg ^= 1;
531 }
532
533 vlong
534 mpgetfix(Mpint *a)
535 {
536 vlong v;
537
538 if(a->ovf) {
539 yyerror("constant overflow");
540 return 0;
541 }
542
543 v = (vlong)a->a[0];
544 v |= (vlong)a->a[1] << Mpscale;
545 v |= (vlong)a->a[2] << (Mpscale+Mpscale);
546 if(a->neg)
547 v = -v;
548 return v;
549 }
550
551 void
552 mpmovecfix(Mpint *a, vlong c)
553 {
554 int i;
555 long *a1;
556 vlong x;
557
558 a->neg = 0;
559 a->ovf = 0;
560
561 x = c;
562 if(x < 0) {
563 a->neg = 1;
564 x = -x;
565 }
566
567 a1 = &a->a[0];
568 for(i=0; i<Mpprec; i++) {
569 *a1++ = x&Mpmask;
570 x >>= Mpscale;
571 }
572 }
573
574 void
575 mpdivmodfixfix(Mpint *q, Mpint *r, Mpint *n, Mpint *d)
576 {
577 int i, ns, ds;
578
579 ns = n->neg;
580 ds = d->neg;
581 n->neg = 0;
582 d->neg = 0;
583
584 mpmovefixfix(r, n);
585 mpmovecfix(q, 0);
586
587 // shift denominator until it
588 // is larger than numerator
589 for(i=0; i<Mpprec*Mpscale; i++) {
590 if(mpcmp(d, r) > 0)
591 break;
592 mplsh(d);
593 }
594
595 // if it never happens
596 // denominator is probably zero
597 if(i >= Mpprec*Mpscale) {
598 q->ovf = 1;
599 r->ovf = 1;
600 n->neg = ns;
601 d->neg = ds;
602 yyerror("set ovf in mpdivmodfixfix");
603 return;
604 }
605
606 // shift denominator back creating
607 // quotient a bit at a time
608 // when done the remaining numerator
609 // will be the remainder
610 for(; i>0; i--) {
611 mplsh(q);
612 mprsh(d);
613 if(mpcmp(d, r) <= 0) {
614 mpaddcfix(q, 1);
615 mpsubfixfix(r, d);
616 }
617 }
618
619 n->neg = ns;
620 d->neg = ds;
621 r->neg = ns;
622 q->neg = ns^ds;
623 }
624
625 static int
626 iszero(Mpint *a)
627 {
628 long *a1;
629 int i;
630 a1 = &a->a[0] + Mpprec;
631 for(i=0; i<Mpprec; i++) {
632 if(*--a1 != 0)
633 return 0;
634 }
635 return 1;
636 }
637
638 void
639 mpdivfract(Mpint *a, Mpint *b)
640 {
641 Mpint n, d;
642 int i, j, neg;
643 long *a1, x;
644
645 mpmovefixfix(&n, a); // numerator
646 mpmovefixfix(&d, b); // denominator
647 a1 = &a->a[Mpprec]; // quotient
648
649 neg = n.neg ^ d.neg;
650 n.neg = 0;
651 d.neg = 0;
652 for(i=0; i<Mpprec; i++) {
653 x = 0;
654 for(j=0; j<Mpscale; j++) {
655 x <<= 1;
656 if(mpcmp(&d, &n) <= 0) {
657 if(!iszero(&d))
658 x |= 1;
659 mpsubfixfix(&n, &d);
660 }
661 mprsh(&d);
662 }
663 *--a1 = x;
664 }
665 a->neg = neg;
666 }
667
668 int
669 mptestfix(Mpint *a)
670 {
671 Mpint b;
672 int r;
673
674 mpmovecfix(&b, 0);
675 r = mpcmp(a, &b);
676 if(a->neg) {
677 if(r > 0)
678 return -1;
679 if(r < 0)
680 return +1;
681 }
682 return r;
683 }