Download

Reference Manual

Navigation

← Previous Changeset
Next Changeset →

Changeset 3884

Timestamp:

08/18/08 03:09:56 (4 months ago)

Author:

Don Clugston

Message:

* multibyteShl() now returns the msb.
* bugfix in Karatsuba multiply for an asymmetric case exactly at the Karatsuba threshold.

Files:

trunk/tango/math/internal/BignumNoAsm.d (modified) (2 diffs)
trunk/tango/math/internal/BignumX86.d (modified) (2 diffs)
trunk/tango/math/internal/BiguintCore.d (modified) (8 diffs)

Legend:

: Unmodified
: Added
: Removed
: Modified
: Copied
: Moved

trunk/tango/math/internal/BignumNoAsm.d

r3831	r3884
156	156	* numbits must be in the range 1..31
157	157	*/
158		~~void~~ multibyteShl(uint [] dest, uint [] src, uint numbits)
	158	uint multibyteShl(uint [] dest, uint [] src, uint numbits)
159	159	{
160	160	ulong c = 0;
…	…
164	164	c >>>= 32;
165	165	}
	166	return cast(uint)c;
166	167	}
167	168

trunk/tango/math/internal/BignumX86.d

r3877	r3884
298	298	/** dest[#] = src[#] << numbits
299	299	* numbits must be in the range 1..31
300		*/
301		void multibyteShl(uint [] dest, uint [] src, uint numbits)
	300	* Returns the overflow
	301	*/
	302	uint multibyteShl(uint [] dest, uint [] src, uint numbits)
302	303	{
303	304	// Timing: Optimal for P6 family.
…	…
333	334	jg L_even;
334	335	L_last:
335		shl EAX, CL;
	336	mov EDX, 0;
	337	shld EAX, EDX, CL;
336	338	mov [EDI], EAX;
337
	339	mov EAX, EDX;
338	340	pop EBX;
339	341	pop EDI;

trunk/tango/math/internal/BiguintCore.d

r3877	r3884
19	19	}
20	20
	21	private:
	22	uint [1] BIGINTZERO = [0];
21	23	public:
22		// Converts a big uint to a hexadecimal string.
23		//
24		// Optionally, a separator character (eg, an underscore) may be added between
25		// every 8 digits.
26		// buff.length must be data.length*8 if separator is zero,
27		// or data.length*9 if separator is non-zero. It will be completely filled.
28		char [] biguintToHex(char [] buff, uint [] data, char separator=0)
29		{
30		int x=0;
31		for (int i=data.length - 1; i>=0; --i) {
32		toHexZeroPadded(buff[x..x+8], data[i]);
33		x+=8;
34		if (separator) {
35		if (i>0) buff[x]='_';
36		++x;
37		}
38		}
39		return buff;
40		}
41
42		/** Convert a big uint into a decimal string.
43		*
44		* Params:
45		* data The biguint to be converted. Will be destroyed.
46		* buff The destination buffer for the decimal string. Must be
47		* large enough to store the result, including leading zeros.
48		* Will be filled starting from the end.
49		*
50		* Ie, buff.length must be (data.length32)/log2(10) = 9.63296 data.length.
51		* Returns:
52		* the lowest index of buff which was used.
53		*/
54		int biguintToDecimal(char [] buff, uint [] data){
55		int sofar = buff.length;
56		// Might be better to divide by (10^38/2^32) since that gives 38 digits for
57		// the price of 3 divisions and a shr; this version only gives 27 digits
58		// for 3 divisions.
59		while(data.length>1) {
60		uint rem = multibyteDivAssign(data, 10_0000_0000, 0);
61		itoaZeroPadded(buff[sofar-9 .. sofar], rem);
62		sofar -= 9;
63		if (data[$-1]==0 && data.length>1) {
64		data.length = data.length - 1;
65		}
66		}
67		itoaZeroPadded(buff[sofar-10 .. sofar], data[0]);
68		sofar-=10;
69		// and strip off the leading zeros
70		while(sofar!= buff.length-1 && buff[sofar] == '0') sofar++;
71		return sofar;
72		}
73
74		/** Convert a decimal string into a big uint.
75		*
76		* Params:
77		* data The biguint to be receive the result. Must be large enough to
78		* store the result.
79		* s The decimal string. May contain 0..9, or _. Will be preserved.
80		*
81		* The required length for the destination buffer is slightly less than
82		* 1 + s.length/log2(10) = 1 + s.length/3.3219.
83		*
84		* Returns:
85		* the highest index of data which was used.
86		*/
87		int biguintFromDecimal(uint [] data, char [] s) {
88		// Convert to base 1e19 = 10_000_000_000_000_000_000.
89		// (this is the largest power of 10 that will fit into a long).
90		// The length will be less than 1 + s.length/log2(10) = 1 + s.length/3.3219.
91		// 485 bits will only just fit into 146 decimal digits.
92		uint lo = 0;
93		uint x = 0;
94		ulong y = 0;
95		uint hi = 0;
96		data[0] = 0; // initially number is 0.
97
98		for (int i= (s[0]=='-' \|\| s[0]=='+')? 1 : 0; i<s.length; ++i) {
99		if (s[i] == '_') continue;
100		x *= 10;
101		x += s[i] - '0';
102		++lo;
103		if (lo==9) {
104		y = x;
105		x = 0;
106		}
107		if (lo==18) {
108		y *= 10_0000_0000;
109		y += x;
110		x = 0;
111		}
112		if (lo==19) {
113		y *= 10;
114		y += x;
115		x = 0;
116		// Multiply existing number by 10^19, then add y1.
117		if (hi>0) {
118		data[hi] = multibyteMul(data[0..hi], data[0..hi], 12207031252, 0); // 5^132 = 0x9184_E72A
119		++hi;
120		data[hi] = multibyteMul(data[0..hi], data[0..hi], 15625262144, 0); // 5^62^18 = 0xF424_0000
121		++hi;
122		} else hi = 2;
123		uint c = multibyteIncrementAssign!('+')(data[0..hi], cast(uint)(y&0xFFFF_FFFF));
124		c += multibyteIncrementAssign!('+')(data[1..hi], cast(uint)(y>>32));
125		if (c!=0) { data[hi]=c; ++hi; }
126		y = 0;
127		lo = 0;
128		}
129		}
130		// Now set y = all remaining digits.
131		if (lo>=18) {
132		} else if (lo>=9) {
133		for (int k=9; k<lo; ++k) y*=10;
134		y+=x;
	24
	25	// return 1 if x>y, -1 if x<y, 0 if equal
	26	int biguintCompare(uint [] x, uint []y)
	27	{
	28	uint k = lastDifferentDigit(x, y);
	29	if (x[k] == y[k]) return 0;
	30	return x[k] > y[k] ? 1 : -1;
	31	}
	32
	33	// return x >> y
	34	uint [] biguintShr(uint[] x, ulong y)
	35	{
	36	assert(y>0);
	37	uint bits = cast(uint)y & 31;
	38	if ((y>>5) >= x.length) return BIGINTZERO;
	39	uint words = cast(uint)(y >> 5);
	40	if (bits==0) {
	41	return x[words..$];
135	42	} else {
136		for (int k=0; k<lo; ++k) y*=10;
137		y+=x;
138		}
139		if (y!=0) {
140		if (hi==0) { cast(ulong )(&data[hi]) = y; hi=2; }
141		else {
142
143		while (lo>0) {
144		uint c = multibyteMul(data[0..hi], data[0..hi], 10, 0);
145		if (c!=0) { data[hi]=c; ++hi; }
146		--lo;
147		}
148		uint c = multibyteIncrementAssign!('+')(data[0..hi], cast(uint)(y&0xFFFF_FFFF));
149		if (y>0xFFFF_FFFFL) {
150		c += multibyteIncrementAssign!('+')(data[1..hi], cast(uint)(y>>32));
151		}
152		if (c!=0) { data[hi]=c; ++hi; }
153		hi+=2;
154		}
155		}
156		return hi;
157		}
158
159		public:
	43	uint [] result = new uint[x.length - words];
	44	multibyteShr(result, x[words..$], bits);
	45	if (result.length>1 && result[$-1]==0) return result[0..$-1];
	46	else return result;
	47	}
	48	}
	49
	50	// return x << y
	51	uint [] biguintShl(uint[] x, ulong y)
	52	{
	53	assert(y>0);
	54	if (x.length==1 && x[0]==0) return x;
	55	uint bits = cast(uint)y & 31;
	56	assert ((y>>5) < cast(ulong)(uint.max));
	57	uint words = cast(uint)(y >> 5);
	58	uint [] result = new uint[x.length + words+1];
	59	result[0..words] = 0;
	60	if (bits==0) {
	61	result[words..words+x.length] = x[];
	62	return result[0..words+x.length];
	63	} else {
	64	uint c = multibyteShl(result[words..words+x.length], x, bits);
	65	if (c==0) return result[0..words+x.length];
	66	result[$-1] = c;
	67	return result;
	68	}
	69	}
	70
160	71
161	72	/** General unsigned subtraction routine for bigints.
…	…
198	109	}
199	110
200		// return a+b
	111	// return a + b
201	112	uint [] biguintAdd(uint[] a, uint [] b) {
202	113	uint [] x, y;
…	…
399	310	}
400	311
	312	public:
	313	// Converts a big uint to a hexadecimal string.
	314	//
	315	// Optionally, a separator character (eg, an underscore) may be added between
	316	// every 8 digits.
	317	// buff.length must be data.length*8 if separator is zero,
	318	// or data.length*9 if separator is non-zero. It will be completely filled.
	319	char [] biguintToHex(char [] buff, uint [] data, char separator=0)
	320	{
	321	int x=0;
	322	for (int i=data.length - 1; i>=0; --i) {
	323	toHexZeroPadded(buff[x..x+8], data[i]);
	324	x+=8;
	325	if (separator) {
	326	if (i>0) buff[x]='_';
	327	++x;
	328	}
	329	}
	330	return buff;
	331	}
	332
	333	/** Convert a big uint into a decimal string.
	334	*
	335	* Params:
	336	* data The biguint to be converted. Will be destroyed.
	337	* buff The destination buffer for the decimal string. Must be
	338	* large enough to store the result, including leading zeros.
	339	* Will be filled starting from the end.
	340	*
	341	* Ie, buff.length must be (data.length32)/log2(10) = 9.63296 data.length.
	342	* Returns:
	343	* the lowest index of buff which was used.
	344	*/
	345	int biguintToDecimal(char [] buff, uint [] data){
	346	int sofar = buff.length;
	347	// Might be better to divide by (10^38/2^32) since that gives 38 digits for
	348	// the price of 3 divisions and a shr; this version only gives 27 digits
	349	// for 3 divisions.
	350	while(data.length>1) {
	351	uint rem = multibyteDivAssign(data, 10_0000_0000, 0);
	352	itoaZeroPadded(buff[sofar-9 .. sofar], rem);
	353	sofar -= 9;
	354	if (data[$-1]==0 && data.length>1) {
	355	data.length = data.length - 1;
	356	}
	357	}
	358	itoaZeroPadded(buff[sofar-10 .. sofar], data[0]);
	359	sofar-=10;
	360	// and strip off the leading zeros
	361	while(sofar!= buff.length-1 && buff[sofar] == '0') sofar++;
	362	return sofar;
	363	}
	364
	365	/** Convert a decimal string into a big uint.
	366	*
	367	* Params:
	368	* data The biguint to be receive the result. Must be large enough to
	369	* store the result.
	370	* s The decimal string. May contain 0..9, or _. Will be preserved.
	371	*
	372	* The required length for the destination buffer is slightly less than
	373	* 1 + s.length/log2(10) = 1 + s.length/3.3219.
	374	*
	375	* Returns:
	376	* the highest index of data which was used.
	377	*/
	378	int biguintFromDecimal(uint [] data, char [] s) {
	379	// Convert to base 1e19 = 10_000_000_000_000_000_000.
	380	// (this is the largest power of 10 that will fit into a long).
	381	// The length will be less than 1 + s.length/log2(10) = 1 + s.length/3.3219.
	382	// 485 bits will only just fit into 146 decimal digits.
	383	uint lo = 0;
	384	uint x = 0;
	385	ulong y = 0;
	386	uint hi = 0;
	387	data[0] = 0; // initially number is 0.
	388
	389	for (int i= (s[0]=='-' \|\| s[0]=='+')? 1 : 0; i<s.length; ++i) {
	390	if (s[i] == '_') continue;
	391	x *= 10;
	392	x += s[i] - '0';
	393	++lo;
	394	if (lo==9) {
	395	y = x;
	396	x = 0;
	397	}
	398	if (lo==18) {
	399	y *= 10_0000_0000;
	400	y += x;
	401	x = 0;
	402	}
	403	if (lo==19) {
	404	y *= 10;
	405	y += x;
	406	x = 0;
	407	// Multiply existing number by 10^19, then add y1.
	408	if (hi>0) {
	409	data[hi] = multibyteMul(data[0..hi], data[0..hi], 12207031252, 0); // 5^132 = 0x9184_E72A
	410	++hi;
	411	data[hi] = multibyteMul(data[0..hi], data[0..hi], 15625262144, 0); // 5^62^18 = 0xF424_0000
	412	++hi;
	413	} else hi = 2;
	414	uint c = multibyteIncrementAssign!('+')(data[0..hi], cast(uint)(y&0xFFFF_FFFF));
	415	c += multibyteIncrementAssign!('+')(data[1..hi], cast(uint)(y>>32));
	416	if (c!=0) { data[hi]=c; ++hi; }
	417	y = 0;
	418	lo = 0;
	419	}
	420	}
	421	// Now set y = all remaining digits.
	422	if (lo>=18) {
	423	} else if (lo>=9) {
	424	for (int k=9; k<lo; ++k) y*=10;
	425	y+=x;
	426	} else {
	427	for (int k=0; k<lo; ++k) y*=10;
	428	y+=x;
	429	}
	430	if (y!=0) {
	431	if (hi==0) { cast(ulong )(&data[hi]) = y; hi=2; }
	432	else {
	433
	434	while (lo>0) {
	435	uint c = multibyteMul(data[0..hi], data[0..hi], 10, 0);
	436	if (c!=0) { data[hi]=c; ++hi; }
	437	--lo;
	438	}
	439	uint c = multibyteIncrementAssign!('+')(data[0..hi], cast(uint)(y&0xFFFF_FFFF));
	440	if (y>0xFFFF_FFFFL) {
	441	c += multibyteIncrementAssign!('+')(data[1..hi], cast(uint)(y>>32));
	442	}
	443	if (c!=0) { data[hi]=c; ++hi; }
	444	hi+=2;
	445	}
	446	}
	447	return hi;
	448	}
	449
401	450
402	451	private:
…	…
465	514	void simpleAddAssign(uint [] result, uint [] right)
466	515	{
467		assert(result.length > right.length);
	516	assert(result.length >= right.length);
468	517	uint c = multibyteAddSub!('+')(result[0..right.length], result[0..right.length], right, 0);
469	518	if (c) {
	519	assert(result.length > right.length);
470	520	c = multibyteIncrementAssign!('+')(result[right.length .. $], c);
471	521	assert(c==0);
…	…
496	546	assert(result.length == x.length + y.length);
497	547	assert(x.length >= y.length);
498		if (x.length <= KARATSUBALIMIT) {
	548	if (y.length <= KARATSUBALIMIT) {
499	549	return mulSimple(result, x, y);
500	550	}
501	551	// Must be almost square.
502		assert(2 * y.length * y.length >~~= x.length * x.length~~, "Asymmetric Karatsuba");
	552	assert(2 * y.length * y.length > (x.length-1) * (x.length-1), "Asymmetric Karatsuba");
503	553
504	554	// Karatsuba multiply uses the following result:
…	…
528	578	// (Nx1 + x0)(Ny1 + y0) = (NN) x1y1 + x0y0 - N * mid
529	579	// where mid = (x0-x1)*(y0-y1) - x1y1 - x0y0
530		// since then ~~the lengths~~ cannot exceed length N.
	580	// since then mid.length cannot exceed length N.
531	581	uint carry_x = addSimple(xsum, x0, x1);
532	582	uint carry_y = addSimple(ysum, y0, y1);
…	…
536	586	if (carry_x) simpleAddAssign(mid[half..$], ysum);
537	587	if (carry_y) simpleAddAssign(mid[half..$], xsum);
538
539	588	// Low half of result gets x0 * y0. High half gets x1 * y1
540	589
…	…
555	604	ysmaller?y1:x1[0..quarter], newscratchbuff);
556	605	// Save the part which will be overwritten.
	606	bool ysmaller2 = ((x1.length -quarter) >= y1.length);
557	607	newscratchbuff[0..y1.length] = resultHigh[quarter..quarter+y1.length];
558		mulKaratsuba(resultHigh[quarter..$], ysmaller?x1[quarter..$]: y1,
559		ysmaller?y1:x1[quarter..$], newscratchbuff[y1.length..$]);
	608	mulKaratsuba(resultHigh[quarter..$], ysmaller2?x1[quarter..$]: y1,
	609	ysmaller2?y1:x1[quarter..$], newscratchbuff[y1.length..$]);
	610
560	611	resultHigh[quarter..$].simpleAddAssign(newscratchbuff[0..y1.length]);
561	612	}
562	613	} else mulKaratsuba(resultHigh, x1, y1, newscratchbuff);
563
	614
564	615	mid.simpleSubAssign(resultHigh);
565
	616
566	617	// result += MID
567	618	result[half..$].simpleAddAssign(mid);
568	619	}
	620
569	621
570	622	private:

Download in other formats:

Trac

dsource.org

Copyright © 2006-2008 Tango. All Rights Reserved. | Page Width: Static or Dynamic