plan9port

fork of plan9port with libvec, libstr and libsdb
Log | Files | Refs | README | LICENSE

mpaux.c (2667B)


      1 #include "os.h"
      2 #include <mp.h>
      3 #include "dat.h"
      4 
      5 static mpdigit _mptwodata[1] = { 2 };
      6 static mpint _mptwo =
      7 {
      8 	1,
      9 	1,
     10 	1,
     11 	_mptwodata,
     12 	MPstatic
     13 };
     14 mpint *mptwo = &_mptwo;
     15 
     16 static mpdigit _mponedata[1] = { 1 };
     17 static mpint _mpone =
     18 {
     19 	1,
     20 	1,
     21 	1,
     22 	_mponedata,
     23 	MPstatic
     24 };
     25 mpint *mpone = &_mpone;
     26 
     27 static mpdigit _mpzerodata[1] = { 0 };
     28 static mpint _mpzero =
     29 {
     30 	1,
     31 	1,
     32 	0,
     33 	_mpzerodata,
     34 	MPstatic
     35 };
     36 mpint *mpzero = &_mpzero;
     37 
     38 static int mpmindigits = 33;
     39 
     40 /* set minimum digit allocation */
     41 void
     42 mpsetminbits(int n)
     43 {
     44 	if(n < 0)
     45 		sysfatal("mpsetminbits: n < 0");
     46 	if(n == 0)
     47 		n = 1;
     48 	mpmindigits = DIGITS(n);
     49 }
     50 
     51 /* allocate an n bit 0'd number  */
     52 mpint*
     53 mpnew(int n)
     54 {
     55 	mpint *b;
     56 
     57 	if(n < 0)
     58 		sysfatal("mpsetminbits: n < 0");
     59 
     60 	b = mallocz(sizeof(mpint), 1);
     61 	if(b == nil)
     62 		sysfatal("mpnew: %r");
     63 	n = DIGITS(n);
     64 	if(n < mpmindigits)
     65 		n = mpmindigits;
     66 	b->p = (mpdigit*)mallocz(n*Dbytes, 1);
     67 	if(b->p == nil)
     68 		sysfatal("mpnew: %r");
     69 	b->size = n;
     70 	b->sign = 1;
     71 
     72 	return b;
     73 }
     74 
     75 /* guarantee at least n significant bits */
     76 void
     77 mpbits(mpint *b, int m)
     78 {
     79 	int n;
     80 
     81 	n = DIGITS(m);
     82 	if(b->size >= n){
     83 		if(b->top >= n)
     84 			return;
     85 		memset(&b->p[b->top], 0, Dbytes*(n - b->top));
     86 		b->top = n;
     87 		return;
     88 	}
     89 	b->p = (mpdigit*)realloc(b->p, n*Dbytes);
     90 	if(b->p == nil)
     91 		sysfatal("mpbits: %r");
     92 	memset(&b->p[b->top], 0, Dbytes*(n - b->top));
     93 	b->size = n;
     94 	b->top = n;
     95 }
     96 
     97 void
     98 mpfree(mpint *b)
     99 {
    100 	if(b == nil)
    101 		return;
    102 	if(b->flags & MPstatic)
    103 		sysfatal("freeing mp constant");
    104 	memset(b->p, 0, b->size*Dbytes);	/* information hiding */
    105 	free(b->p);
    106 	free(b);
    107 }
    108 
    109 void
    110 mpnorm(mpint *b)
    111 {
    112 	int i;
    113 
    114 	for(i = b->top-1; i >= 0; i--)
    115 		if(b->p[i] != 0)
    116 			break;
    117 	b->top = i+1;
    118 	if(b->top == 0)
    119 		b->sign = 1;
    120 }
    121 
    122 mpint*
    123 mpcopy(mpint *old)
    124 {
    125 	mpint *new;
    126 
    127 	new = mpnew(Dbits*old->size);
    128 	new->top = old->top;
    129 	new->sign = old->sign;
    130 	memmove(new->p, old->p, Dbytes*old->top);
    131 	return new;
    132 }
    133 
    134 void
    135 mpassign(mpint *old, mpint *new)
    136 {
    137 	mpbits(new, Dbits*old->top);
    138 	new->sign = old->sign;
    139 	new->top = old->top;
    140 	memmove(new->p, old->p, Dbytes*old->top);
    141 }
    142 
    143 /* number of significant bits in mantissa */
    144 int
    145 mpsignif(mpint *n)
    146 {
    147 	int i, j;
    148 	mpdigit d;
    149 
    150 	if(n->top == 0)
    151 		return 0;
    152 	for(i = n->top-1; i >= 0; i--){
    153 		d = n->p[i];
    154 		for(j = Dbits-1; j >= 0; j--){
    155 			if(d & (((mpdigit)1)<<j))
    156 				return i*Dbits + j + 1;
    157 		}
    158 	}
    159 	return 0;
    160 }
    161 
    162 /* k, where n = 2**k * q for odd q */
    163 int
    164 mplowbits0(mpint *n)
    165 {
    166 	int k, bit, digit;
    167 	mpdigit d;
    168 
    169 	if(n->top==0)
    170 		return 0;
    171 	k = 0;
    172 	bit = 0;
    173 	digit = 0;
    174 	d = n->p[0];
    175 	for(;;){
    176 		if(d & (1<<bit))
    177 			break;
    178 		k++;
    179 		bit++;
    180 		if(bit==Dbits){
    181 			if(++digit >= n->top)
    182 				return 0;
    183 			d = n->p[digit];
    184 			bit = 0;
    185 		}
    186 	}
    187 	return k;
    188 }