--- a/LargeInteger.st Mon Mar 02 14:28:02 2020 +0100
+++ b/LargeInteger.st Mon Mar 02 14:50:33 2020 +0100
@@ -333,7 +333,57 @@
"Modified: / 8.5.1998 / 21:40:41 / cg"
! !
-
+!LargeInteger class methodsFor:'accessing'!
+
+javaArrayClass
+ ^ SignedLongIntegerArray
+
+ "Created: / 11-02-2011 / 10:51:32 / Jan Vrany <jan.vrany@fit.cvut.cz>"
+!
+
+javaName
+
+ ^'long'.
+
+ "Modified: / 25-02-2011 / 18:59:31 / Jan Vrany <jan.vrany@fit.cvut.cz>"
+!
+
+javaWrapperClass
+ ^(Java classForName: 'java.lang.Long')
+
+ "Created: / 24-02-2012 / 19:42:20 / Jan Vrany <jan.vrany@fit.cvut.cz>"
+! !
+
+!LargeInteger class methodsFor:'autoboxing support'!
+
+javaBox: anObject
+ | wrapper |
+
+ wrapper := (JavaVM classForName: 'java.lang.Long') new.
+ wrapper perform: #'<init>(J)V' with: anObject with: nil.
+ ^ wrapper
+
+ "Created: / 16-08-2011 / 09:58:40 / Jan Vrany <jan.vrany@fit.cvut.cz>"
+ "Modified: / 14-11-2013 / 01:21:21 / Jan Vrany <jan.vrany@fit.cvut.cz>"
+!
+
+javaUnbox: object onError: errorBlock
+
+ | value |
+
+ (object class binaryName = 'java/lang/Long') ifFalse:[
+ errorBlock value.
+ ].
+
+ value := object instVarNamed:#value.
+ (value between: "Integer.MIN_VALUE"-9223372036854775808 and: "Integer.MAX_VALUE" 9223372036854775807) ifFalse:[
+ errorBlock value.
+ ].
+ ^value
+
+ "Created: / 22-11-2011 / 11:45:30 / Jan Vrany <jan.vrany@fit.cvut.cz>"
+ "Modified: / 08-10-2013 / 22:42:10 / Jan Vrany <jan.vrany@fit.cvut.cz>"
+! !
!LargeInteger class methodsFor:'coercing & converting'!
@@ -352,6 +402,13 @@
^ true
"Modified: 23.4.1996 / 15:59:21 / cg"
+!
+
+isJavaPrimitiveType
+
+ ^true
+
+ "Created: / 04-02-2011 / 11:55:50 / Jan Vrany <jan.vrany@fit.cvut.cz>"
! !
!LargeInteger methodsFor:'arithmetic'!
@@ -4715,8 +4772,8 @@
otherDigitByteArray := aLargeInteger digitBytes.
len2 := otherDigitByteArray size.
lenRslt := len1 + len2.
- UseKarazuba ~~ false ifTrue:[
- lenRslt > 400 ifTrue:[ ^ self absMulKarazuba:aLargeInteger ].
+ (UseKarazuba ~~ false and:[lenRslt > 900]) ifTrue:[
+ ^ self absMulKarazuba:aLargeInteger.
].
result := self class basicNew numberOfDigits:lenRslt.
@@ -4726,239 +4783,250 @@
if (__isByteArray(__INST(digitByteArray))
&& __isByteArray(otherDigitByteArray)
&& __isByteArray(resultDigitByteArray)) {
- unsigned char *myBytes = __ByteArrayInstPtr(__INST(digitByteArray))->ba_element;
- unsigned char *otherBytes = __ByteArrayInstPtr(otherDigitByteArray)->ba_element;
- unsigned char *resultBytes = __ByteArrayInstPtr(resultDigitByteArray)->ba_element;
- unsigned char *_p1, *_pResult0, *_p1Last, *_p2Last;
- unsigned INT _v;
- int _len1 = __intVal(len1);
- int _len2 = __intVal(len2);
-
- _p1Last = myBytes + _len1 - 1; /* the last byte */
- _p2Last = otherBytes + _len2 - 1; /* the last byte */
- _pResult0 = resultBytes;
-
- /*
- * aaa...aaa f1[0] * f2
- * + bbb...bbb f1[1] * f2
- * + ccc...ccc f1[2] * f2
- * + ...
- * + xxx...xxx f1[high] * f2
- *
- * start short-wise
- * bounds: (16rFFFF * 16rFFFF) + 16rFFFF -> FFFF0000
- */
- _p1 = myBytes;
+ unsigned char *myBytes = __ByteArrayInstPtr(__INST(digitByteArray))->ba_element;
+ unsigned char *otherBytes = __ByteArrayInstPtr(otherDigitByteArray)->ba_element;
+ unsigned char *resultBytes = __ByteArrayInstPtr(resultDigitByteArray)->ba_element;
+ unsigned char *_p1, *_pResult0, *_p1Last, *_p2Last;
+ unsigned INT _v;
+ int _len1 = __intVal(len1);
+ int _len2 = __intVal(len2);
+
+ _p1Last = myBytes + _len1 - 1; /* the last byte */
+ _p2Last = otherBytes + _len2 - 1; /* the last byte */
+ _pResult0 = resultBytes;
+
+ /*
+ * aaa...aaa f1[0] * f2
+ * + bbb...bbb f1[1] * f2
+ * + ccc...ccc f1[2] * f2
+ * + ...
+ * + xxx...xxx f1[high] * f2
+ *
+ * start short-wise
+ * bounds: (16rFFFF * 16rFFFF) + 16rFFFF -> FFFF0000
+ */
+ _p1 = myBytes;
#if defined(__LSBFIRST__) && (__POINTER_SIZE__ == 8)
- /* loop over ints of f1 */
- for (; _p1 <= (_p1Last-3); _p1 += 4, _pResult0 += 4) {
- unsigned char *_pResult = _pResult0;
- unsigned char *_p2;
- unsigned INT word1 = ((unsigned int *)_p1)[0];
-
- _v = 0;
-
- /* loop over ints of f2 */
- for (_p2 = otherBytes; _p2 <= (_p2Last-3); _p2 += 4) {
- _v = word1 * (unsigned INT)(((unsigned int *)_p2)[0])
- + _v + (unsigned INT)((unsigned int *)_pResult)[0];
- ((unsigned int *)_pResult)[0] = _v /* & 0xFFFFFFFF */;
- _v >>= 32; /* now _v contains the carry*/
- _pResult += 4;
- }
-
- /* possible odd up to 3 highBytes of f2 */
- for ( ; _p2 <= _p2Last; _p2++) {
- _v = word1 * (unsigned INT)(_p2[0])
- + _v + (unsigned INT)(_pResult[0]);
-
- ((unsigned char *)_pResult)[0] = _v /* & 0xFFFFFFFF */;
- _v >>= 8; /* now _v contains the carry*/
- _pResult++;
- }
- /* distribute leftover carry byte-wise */
- for ( ; _v; _v >>= 8, _pResult++) {
- _v += _pResult[0];
- _pResult[0] = _v /* & 0xFF */;
- }
- }
+ /* loop over ints of f1 */
+ for (; _p1 <= (_p1Last-3); _p1 += 4, _pResult0 += 4) {
+ unsigned char *_pResult = _pResult0;
+ unsigned char *_p2;
+ unsigned INT word1 = ((unsigned int *)_p1)[0];
+
+ _v = 0;
+
+ /* loop over ints of f2 */
+ for (_p2 = otherBytes; _p2 <= (_p2Last-3); _p2 += 4) {
+ _v = word1 * (unsigned INT)(((unsigned int *)_p2)[0])
+ + _v + (unsigned INT)((unsigned int *)_pResult)[0];
+ ((unsigned int *)_pResult)[0] = _v /* & 0xFFFFFFFF */;
+ _v >>= 32; /* now _v contains the carry*/
+ _pResult += 4;
+ }
+
+ /* possible odd up to 3 highBytes of f2 */
+ for ( ; _p2 <= _p2Last; _p2++) {
+ _v = word1 * (unsigned INT)(_p2[0])
+ + _v + (unsigned INT)(_pResult[0]);
+
+ ((unsigned char *)_pResult)[0] = _v /* & 0xFFFFFFFF */;
+ _v >>= 8; /* now _v contains the carry*/
+ _pResult++;
+ }
+ /* distribute leftover carry byte-wise */
+ for ( ; _v; _v >>= 8, _pResult++) {
+ _v += _pResult[0];
+ _pResult[0] = _v /* & 0xFF */;
+ }
+ }
#endif /* 64bit */
- /* possible odd high short of f1 (or shortLoop, if not 64bit) */
-
- for (; _p1 < _p1Last; _p1 += 2, _pResult0 += 2) {
- unsigned char *_pResult = _pResult0;
- unsigned char *_p2 = otherBytes;
- unsigned int short1 = ((unsigned short *)_p1)[0];
+ /* possible odd high short of f1 (or shortLoop, if not 64bit) */
+
+ for (; _p1 < _p1Last; _p1 += 2, _pResult0 += 2) {
+ unsigned char *_pResult = _pResult0;
+ unsigned char *_p2 = otherBytes;
+ unsigned int short1 = ((unsigned short *)_p1)[0];
#if !defined(__LSBFIRST__)
- short1 = ((short1 >> 8) & 0xFF) | ((short1 & 0xFF) << 8);
+ short1 = ((short1 >> 8) & 0xFF) | ((short1 & 0xFF) << 8);
#endif
- _v = 0;
-
- /* loop over shorts of f2 */
- for ( ; _p2 < _p2Last; _p2 += 2, _pResult += 2) {
+ _v = 0;
+
+ /* loop over shorts of f2 */
+ for ( ; _p2 < _p2Last; _p2 += 2, _pResult += 2) {
#if !defined(__LSBFIRST__)
- unsigned int _short2 = ((unsigned short *)_p2)[0];
- unsigned int _short3 = ((unsigned short *)_pResult)[0];
-
- _short2 = ((_short2 >> 8) /* & 0xFF */) | ((_short2 & 0xFF) << 8);
- _short3 = ((_short3 >> 8) /* & 0xFF */) | ((_short3 & 0xFF) << 8);
- _v = (short1 * _short2) + _short3 + _v;
- _pResult[0] = _v;
- _pResult[1] = _v >> 8;
+ unsigned int _short2 = ((unsigned short *)_p2)[0];
+ unsigned int _short3 = ((unsigned short *)_pResult)[0];
+
+ _short2 = ((_short2 >> 8) /* & 0xFF */) | ((_short2 & 0xFF) << 8);
+ _short3 = ((_short3 >> 8) /* & 0xFF */) | ((_short3 & 0xFF) << 8);
+ _v = (short1 * _short2) + _short3 + _v;
+ _pResult[0] = _v;
+ _pResult[1] = _v >> 8;
#else /* __LSBFIRST__ */
- _v = (short1 * ((unsigned short *)_p2)[0]) + _v + ((unsigned short *)_pResult)[0];
- ((unsigned short *)_pResult)[0] = _v /* & 0xFFFF */;
+ _v = (short1 * ((unsigned short *)_p2)[0]) + _v + ((unsigned short *)_pResult)[0];
+ ((unsigned short *)_pResult)[0] = _v /* & 0xFFFF */;
#endif /* __LSBFIRST__ */
- _v >>= 16; /* now _v contains the carry*/
- }
-
- /* possible odd highByte of f2 */
- for ( ; _p2 <= _p2Last; _p2++, _pResult += 2) {
+ _v >>= 16; /* now _v contains the carry*/
+ }
+
+ /* possible odd highByte of f2 */
+ for ( ; _p2 <= _p2Last; _p2++, _pResult += 2) {
#if !defined(__LSBFIRST__)
- unsigned int _short3 = ((unsigned short *)_pResult)[0];
-
- _short3 = ((_short3 >> 8) /* & 0xFF */) | ((_short3 & 0xFF) << 8);
- _v = (short1 * _p2[0]) + _v + _short3;
- _pResult[0] = _v;
- _pResult[1] = _v >> 8;
+ unsigned int _short3 = ((unsigned short *)_pResult)[0];
+
+ _short3 = ((_short3 >> 8) /* & 0xFF */) | ((_short3 & 0xFF) << 8);
+ _v = (short1 * _p2[0]) + _v + _short3;
+ _pResult[0] = _v;
+ _pResult[1] = _v >> 8;
#else /* __LSBFIRST__ */
- _v = (short1 * _p2[0]) + _v + ((unsigned short *)_pResult)[0];
- ((unsigned short *)_pResult)[0] = _v /* & 0xFFFF */;
+ _v = (short1 * _p2[0]) + _v + ((unsigned short *)_pResult)[0];
+ ((unsigned short *)_pResult)[0] = _v /* & 0xFFFF */;
#endif /* __LSBFIRST__ */
- _v >>= 16; /* now _v contains the carry*/
- }
- /* distribute leftover carry byte-wise */
- for ( ; _v; _v >>= 8, _pResult++) {
- _v += _pResult[0];
- _pResult[0] = _v /* & 0xFF */;
- }
- }
-
- /* possible odd highByte of f1 (or byteLoop, if above is ever disabled) */
- for (; _p1 <= _p1Last; _p1++, _pResult0++) {
- unsigned char *_pResult = _pResult0;
- unsigned char *_p2 = otherBytes;
- unsigned int byte1 = _p1[0];
-
- _v = 0;
-
- /* loop over shorts of f2 */
- for ( ; _p2 < _p2Last; _p2 += 2, _pResult += 2) {
+ _v >>= 16; /* now _v contains the carry*/
+ }
+ /* distribute leftover carry byte-wise */
+ for ( ; _v; _v >>= 8, _pResult++) {
+ _v += _pResult[0];
+ _pResult[0] = _v /* & 0xFF */;
+ }
+ }
+
+ /* possible odd highByte of f1 (or byteLoop, if above is ever disabled) */
+ for (; _p1 <= _p1Last; _p1++, _pResult0++) {
+ unsigned char *_pResult = _pResult0;
+ unsigned char *_p2 = otherBytes;
+ unsigned int byte1 = _p1[0];
+
+ _v = 0;
+
+ /* loop over shorts of f2 */
+ for ( ; _p2 < _p2Last; _p2 += 2, _pResult += 2) {
#if !defined(__LSBFIRST__)
- unsigned int _short2 = ((unsigned short *)_p2)[0];
- unsigned int _short3 = ((unsigned short *)_pResult)[0];
-
- _short2 = ((_short2 >> 8) /* & 0xFF */) | ((_short2 & 0xFF) << 8);
- _short3 = ((_short3 >> 8) /* & 0xFF */) | ((_short3 & 0xFF) << 8);
- _v = (byte1 * _short2) + _v +_short3;
- _pResult[0] = _v;
- _pResult[1] = _v >> 8;
+ unsigned int _short2 = ((unsigned short *)_p2)[0];
+ unsigned int _short3 = ((unsigned short *)_pResult)[0];
+
+ _short2 = ((_short2 >> 8) /* & 0xFF */) | ((_short2 & 0xFF) << 8);
+ _short3 = ((_short3 >> 8) /* & 0xFF */) | ((_short3 & 0xFF) << 8);
+ _v = (byte1 * _short2) + _v +_short3;
+ _pResult[0] = _v;
+ _pResult[1] = _v >> 8;
#else /* __LSBFIRST__ */
- _v = (byte1 * ((unsigned short *)_p2)[0]) + _v + ((unsigned short *)_pResult)[0];
- ((unsigned short *)_pResult)[0] = _v /* & 0xFFFF */;
+ _v = (byte1 * ((unsigned short *)_p2)[0]) + _v + ((unsigned short *)_pResult)[0];
+ ((unsigned short *)_pResult)[0] = _v /* & 0xFFFF */;
#endif /* __LSBFIRST__ */
- _v >>= 16; /* now _v contains the carry*/
- }
-
- /* possible odd highByte of f2 (or byteLoop, if not __LSBFIRST__) */
- for ( ; _p2 <= _p2Last; _p2++, _pResult++) {
- _v = (byte1 * _p2[0]) + _v + _pResult[0];
- _pResult[0] = _v /* & 0xFF */;
- _v >>= 8; /* now _v contains the carry*/
- }
- /* distribute leftover carry byte-wise */
- for ( ; _v; _v >>= 8, _pResult++) {
- _v += _pResult[0];
- _pResult[0] = _v /* & 0xFF */;
- }
- }
- ok = true;
+ _v >>= 16; /* now _v contains the carry*/
+ }
+
+ /* possible odd highByte of f2 (or byteLoop, if not __LSBFIRST__) */
+ for ( ; _p2 <= _p2Last; _p2++, _pResult++) {
+ _v = (byte1 * _p2[0]) + _v + _pResult[0];
+ _pResult[0] = _v /* & 0xFF */;
+ _v >>= 8; /* now _v contains the carry*/
+ }
+ /* distribute leftover carry byte-wise */
+ for ( ; _v; _v >>= 8, _pResult++) {
+ _v += _pResult[0];
+ _pResult[0] = _v /* & 0xFF */;
+ }
+ }
+ ok = true;
}
%}.
ok ifFalse:[
- 1 to:len1 do:[:index1 |
- 1 to:len2 do:[:index2 |
- dstIndex := index1 + index2 - 1.
- prod := (digitByteArray basicAt:index1) * (otherDigitByteArray basicAt:index2).
- prod := prod + (resultDigitByteArray basicAt:dstIndex).
- resultDigitByteArray basicAt:dstIndex put:(prod bitAnd:16rFF).
- carry := prod bitShift:-8.
- carry ~~ 0 ifTrue:[
- idx := dstIndex + 1.
- [carry ~~ 0] whileTrue:[
- v := (resultDigitByteArray basicAt:idx) + carry.
- resultDigitByteArray basicAt:idx put:(v bitAnd:255).
- carry := v bitShift:-8.
- idx := idx + 1
- ]
- ]
- ]
- ].
+ 1 to:len1 do:[:index1 |
+ 1 to:len2 do:[:index2 |
+ dstIndex := index1 + index2 - 1.
+ prod := (digitByteArray basicAt:index1) * (otherDigitByteArray basicAt:index2).
+ prod := prod + (resultDigitByteArray basicAt:dstIndex).
+ resultDigitByteArray basicAt:dstIndex put:(prod bitAnd:16rFF).
+ carry := prod bitShift:-8.
+ carry ~~ 0 ifTrue:[
+ idx := dstIndex + 1.
+ [carry ~~ 0] whileTrue:[
+ v := (resultDigitByteArray basicAt:idx) + carry.
+ resultDigitByteArray basicAt:idx put:(v bitAnd:255).
+ carry := v bitShift:-8.
+ idx := idx + 1
+ ]
+ ]
+ ]
+ ].
].
^ result compressed
+
+ "Modified: / 02-03-2020 / 14:42:11 / Stefan Vogel"
!
-absMulKarazuba:aLargeInteger
+absMulKarazuba:anInteger
"return a LargeInteger representing abs(self) * abs(theArgument) using the karazuba algorithm.
- a = (2^n * p) + q
- b = (2^n * r) + s
- a * b = ((2^n * p) + q) * ((2^n * r) + s)
- = 2^(n+n)*p*r + 2^n*p*s + 2^n*q*r + q*s
- = 2^(n+n)*p*r + (p*r + q*s - (q-p)*(s-r))*2^n + q*s
-
- this is faster for sufficient large n1,n2
- because regular multiplication is O(n1*n2) and karazuma multiplies much smaller numbers
- (half number of bits) but does more multiplications (on smaller numbers) and req's more
- additions and subtractions (on smaller numbers).
- The break-even for when to use regular multiplication has been determined heuristically
- and is somewhere around 1600 bits (digitLength of 200).
- (see test in absMul:)
-
- To disable karazuba, set UseKarazuba to false.
+ a = (2^n * p) + q
+ b = (2^n * r) + s
+ a * b = ((2^n * p) + q) * ((2^n * r) + s)
+ = 2^(n+n)*p*r + 2^n*p*s + 2^n*q*r + q*s
+ = 2^(n+n)*p*r + (p*r + q*s - (q-p)*(s-r))*2^n + q*s
+
+ this is faster for sufficient large n1,n2
+ because regular multiplication is O(n1*n2) and karazuma multiplies much smaller numbers
+ (half number of bits) but does more multiplications (on smaller numbers) and req's more
+ additions and subtractions (on smaller numbers).
+ The break-even for when to use regular multiplication has been determined heuristically
+ and is somewhere around 1600 bits (digitLength of 200 for 32-bit (900 for 64-bit).
+ (see test in absMul:)
+
+ To disable karazuba, set UseKarazuba to false.
"
"/ compute half-sized pi and qi...
- |l1 l2 n nh p q r s pr qs q_p s_r q_ps_r pr_raisedTo_n mdl_raisedTo_nh|
-
- l1 := self digitLength.
- l2 := aLargeInteger digitLength.
-
- n := l1 max:l2.
-
- nh := (n / 2) ceiling. n := nh * 2.
- p := LargeInteger digitBytes:(digitByteArray copyFrom:nh+1) sign:1. "/ high bits of nr1
- q := LargeInteger digitBytes:(digitByteArray copyToMax:nh) sign:1. "/ low bits of nr1
- r := LargeInteger digitBytes:(aLargeInteger digitBytes copyFrom:nh+1) sign:1. "/ high bits of nr2
- s := LargeInteger digitBytes:(aLargeInteger digitBytes copyToMax:nh). "/ low bits of nr2
-
- p := p compressed.
- q := q compressed.
- r := r compressed.
- s := s compressed.
+ |len1 "{ Class:SmallInteger }" len2 "{ Class:SmallInteger }"
+ n "{ Class:SmallInteger }" nh "{ Class:SmallInteger }"
+ otherDigitByteArray p q r s pr qs pr_raisedTo_n mdl_raisedTo_nh|
+
+ len1 := digitByteArray size.
+ anInteger class == SmallInteger ifTrue:[
+ anInteger == 0 ifTrue:[
+ ^ 0.
+ ].
+ r := 0. "/ high bits of nr2
+ s := anInteger. "/ low bits of nr2
+ n := len1.
+ nh := len1 // 2.
+ ] ifFalse:[
+ otherDigitByteArray := anInteger digitBytes.
+ len2 := otherDigitByteArray size.
+ n := len1 max:len2.
+ nh := (n / 2) ceiling.
+ n := nh * 2.
+ r := (LargeInteger digitBytes:(otherDigitByteArray copyFrom:nh+1)) compressed. "/ high bits of nr2
+ s := (LargeInteger digitBytes:(otherDigitByteArray copyToMax:nh)) compressed. "/ low bits of nr2
+ ].
+
+ p := (LargeInteger digitBytes:(digitByteArray copyFrom:nh+1)) compressed. "/ high bits of nr1
+ q := (LargeInteger digitBytes:(digitByteArray copyToMax:nh)) compressed. "/ low bits of nr1
pr := p * r.
qs := q * s.
- q_p := q - p.
- s_r := s - r.
- q_ps_r := q_p * s_r.
-
- pr_raisedTo_n := (pr bitShift:(n*8)).
- mdl_raisedTo_nh := ((pr + qs - q_ps_r) bitShift:(nh*8)).
+
+ pr_raisedTo_n := pr bitShift:(n*8).
+ mdl_raisedTo_nh := (pr + qs - ((q-p)*(s-r))) bitShift:(nh*8).
^ pr_raisedTo_n + mdl_raisedTo_nh + qs
"
#(100 500 1000 2500 5000 10000 25000 50000 100000 250000 500000 1000000) do:[:exp |
- |nr r1 r2|
- nr := (3 raisedTo:exp) asInteger.
- Transcript show:exp; show:' nbytes: '; showCR:nr digitLength;
- show:' normal: '; show:(Time microsecondsToRun:[ UseKarazuba := false. r1 := nr * nr ]); showCR:'us';
- show:' karazuba: '; show:(Time microsecondsToRun:[ UseKarazuba := true. r2 := nr absMulKarazuba: nr]); showCR:'us'.
- self assert:(r1 = r2).
+ |nr r1 r2|
+ nr := (3 raisedTo:exp) + 1111111.
+ Transcript show:exp; show:' nbytes: '; showCR:nr digitLength;
+ show:' normal: '; showCR:(TimeDuration toRun:[ UseKarazuba := false. r1 := nr * nr ]);
+ show:' karazuba: '; showCR:(TimeDuration toRun:[ UseKarazuba := true. r2 := nr absMulKarazuba: nr]).
+ self assert:(r1 = r2).
]
+
+ self assert:((3 raisedTo:10000) absMulKarazuba:11111) = (11111 * (3 raisedTo:10000))
"
+
+ "Modified (comment): / 02-03-2020 / 14:50:22 / Stefan Vogel"
!
absPlus:aLargeInteger sign:newSign