Long multiplication
You are encouraged to solve this task according to the task description, using any language you may know.
In this task, explicitly implement long multiplication. This is one possible approach to arbitrary-precision integer algebra.
For output, display the result of 2^64 * 2^64. The decimal representation of 2^64 is:
18446744073709551616
The output of 2^64 * 2^64 is 2^128, and that is:
340282366920938463463374607431768211456
Ada
The following implementation uses representation of a long number by an array of 32-bit elements: <lang ada>type Long_Number is array (Natural range <>) of Unsigned_32;
function "*" (Left, Right : Long_Number) return Long_Number is
Result : Long_Number (0..Left'Length + Right'Length - 1) := (others => 0); Accum : Unsigned_64;
begin
for I in Left'Range loop
for J in Right'Range loop
Accum := Unsigned_64 (Left (I)) * Unsigned_64 (Right (J));
for K in I + J..Result'Last loop
exit when Accum = 0;
Accum := Accum + Unsigned_64 (Result (K));
Result (K) := Unsigned_32 (Accum and 16#FFFF_FFFF#);
Accum := Accum / 2**32;
end loop;
end loop;
end loop;
for Index in reverse Result'Range loop -- Normalization
if Result (Index) /= 0 then
return Result (0..Index);
end if;
end loop;
return (0 => 0);
end "*";</lang>
The task requires conversion into decimal base. For this we also need division to short number with a remainder. Here it is: <lang ada>procedure Div
( Dividend : in out Long_Number;
Last : in out Natural;
Remainder : out Unsigned_32;
Divisor : Unsigned_32
) is
Div : constant Unsigned_64 := Unsigned_64 (Divisor);
Accum : Unsigned_64 := 0;
Size : Natural := 0;
begin
for Index in reverse Dividend'First..Last loop
Accum := Accum * 2**32 + Unsigned_64 (Dividend (Index));
Dividend (Index) := Unsigned_32 (Accum / Div);
if Size = 0 and then Dividend (Index) /= 0 then
Size := Index;
end if;
Accum := Accum mod Div;
end loop;
Remainder := Unsigned_32 (Accum);
Last := Size;
end Div;</lang>
With the above the test program: <lang ada>with Ada.Strings.Unbounded; use Ada.Strings.Unbounded; with Ada.Text_IO; use Ada.Text_IO; with Interfaces; use Interfaces;
procedure Long_Multiplication is
-- Insert definitions above here
procedure Put (Value : Long_Number) is
X : Long_Number := Value;
Last : Natural := X'Last;
Digit : Unsigned_32;
Result : Unbounded_String;
begin
loop
Div (X, Last, Digit, 10);
Append (Result, Character'Val (Digit + Character'Pos ('0')));
exit when Last = 0 and then X (0) = 0;
end loop;
for Index in reverse 1..Length (Result) loop
Put (Element (Result, Index));
end loop;
end Put;
X : Long_Number := (0 => 0, 1 => 0, 2 => 1) * (0 => 0, 1 => 0, 2 => 1);
begin
Put (X);
end Long_Multiplication;</lang>
Sample output:
340282366920938463463374607431768211456
ALGOL 68
The long multiplication for the golden ratio has been included as half the digits cancel and end up as being zero. This is useful for testing.
Built in or standard distribution routines
ALGOL 68G allows any precision for long long int to be defined when the program is run, e.g. 200 digits. <lang algol68>PRAGMAT precision=200 PRAGMAT MODE INTEGER = LONG LONG INT;
LONG INT default integer width := 69; INT width = 69+2;
INT fix w = 1, fix h = 1; # round up #
LONG LONG INT golden ratio w := ENTIER ((long long sqrt(5)-1) / 2 * LENG LENG 10 ** default integer width + fix w),
golden ratio h := ENTIER ((long long sqrt(5)+1) / 2 * LENG LENG 10 ** default integer width + fix h);
test: (
print(( "The approximate golden ratios, width: ", whole(golden ratio w,width), new line, " length: ", whole(golden ratio h,width), new line, " product is exactly: ", whole(golden ratio w*golden ratio h,width*2), new line));
INTEGER two to the power of 64 = LONG 2 ** 64;
INTEGER neg two to the power of 64 = -(LONG 2 ** 64);
print(("2 ** 64 * -(2 ** 64) = ", whole(two to the power of 64*neg two to the power of 64,width), new line))
)</lang> Output:
The approximate golden ratios, width: +618033988749894848204586834365638117720309179805762862135448622705261
length: +1618033988749894848204586834365638117720309179805762862135448622705261
product is exactly: +1000000000000000000000000000000000000000000000000000000000000000000001201173450350400438606015942314498798603569682901026716145698077078121
2 ** 64 * -(2 ** 64) = -340282366920938463463374607431768211456
Implementation example
<lang algol68>MODE DIGIT = INT; MODE INTEGER = FLEX[0]DIGIT; # an arbitary number of digits #
- "digits" are stored in digit base ten, but 10000 & 2**n (inc hex) can be used #
INT digit base = 1000;
- if possible, then print the digit with one character #
STRING hex digit repr = "0123456789abcdefghijklmnopqrstuvwxyz"[AT 0]; INT digit base digit width = ( digit base <= UPB hex digit repr + 1 | 1 | 1 + ENTIER log(digit base-1) );
INT next digit = -1; # reverse order so digits appear in "normal" order when printed #
PROC raise value error = ([]STRING args)VOID:
( print(("Value Error: ", args, new line)); stop );
PROC raise not implemented error = ([]STRING args)VOID:
( print(("Not implemented Error: ", args, new line)); stop );
PROC raise integer not implemented error = (STRING message)INTEGER:
( raise not implemented error(("INTEGER ", message)); SKIP );
INT half max int = max int OVER 2; IF digit base > half max int THEN raise value error("INTEGER addition may fail") FI;
INT sqrt max int = ENTIER sqrt(max int); IF digit base > sqrt max int THEN raise value error("INTEGER multiplication may fail") FI;
- initialise/cast a INTEGER from a LONG LONG INT #
OP INTEGERINIT = (LONG LONG INT number)INTEGER:(
[1 + ENTIER (SHORTEN SHORTEN long long log(ABS number) / log(digit base))]DIGIT out; LONG LONG INT carry := number; FOR digit out FROM UPB out BY next digit TO LWB out DO LONG LONG INT prev carry := carry; carry %:= digit base; # avoid MOD as it doesn't under handle -ve numbers # out[digit out] := SHORTEN SHORTEN (prev carry - carry * digit base) OD; out
);
- initialise/cast a INTEGER from an LONG INT #
OP INTEGERINIT = (LONG INT number)INTEGER: INTEGERINIT LENG number;
- initialise/cast a INTEGER from an INT #
OP INTEGERINIT = (INT number)INTEGER: INTEGERINIT LENG LENG number;
- remove leading zero "digits" #
OP NORMALISE = ([]DIGIT number)INTEGER: (
INT leading zeros := LWB number - 1; FOR digit number FROM LWB number TO UPB number WHILE number[digit number] = 0 DO leading zeros := digit number OD; IF leading zeros = UPB number THEN 0 ELSE number[leading zeros+1:] FI
);
Define a standard representation for the INTEGER mode. Note: this is rather crude because for a large "digit base" the number is represented as blocks of decimals. It works nicely for powers of ten (10,100,1000,...), but for most larger bases (greater then 35) the repr will be a surprise.
OP REPR = (DIGIT d)STRING:
IF digit base > UPB hex digit repr THEN
STRING out := whole(ABS d, -digit base digit width);
- Replace spaces with zeros #
FOR digit out FROM LWB out TO UPB out DO
IF out[digit out] = " " THEN out[digit out] := "0" FI
OD;
out
ELSE # small enough to represent as ASCII (hex) characters #
hex digit repr[ABS d]
FI;
OP REPR = (INTEGER number)STRING:(
STRING sep = ( digit base digit width > 1 | "," | "" );
INT width := digit base digit width + UPB sep;
[width * UPB number - UPB sep]CHAR out;
INT leading zeros := LWB out - 1;
FOR digit TO UPB number DO
INT start := digit * width - width + 1;
out[start:start+digit base digit width-1] := REPR number[digit];
IF digit base digit width /= 1 & digit /= UPB number THEN
out[start+digit base digit width] := ","
FI
OD;
- eliminate leading zeros #
FOR digit out FROM LWB out TO UPB out WHILE out[digit out] = "0" OR out[digit out] = sep DO leading zeros := digit out OD;
CHAR sign = ( number[1]<0 | "-" | "+" );
- finally return the semi-normalised result #
IF leading zeros = UPB out THEN "0" ELSE sign + out[leading zeros+1:] FI
);</lang>
<lang algol68>################################################################
- Finally Define the required INTEGER multiplication OPerator. #
OP * = (INTEGER a, b)INTEGER:(
- initialise out to all zeros #
[UPB a + UPB b]INT ab; FOR place ab TO UPB ab DO ab[place ab]:=0 OD;
FOR place a FROM UPB a BY next digit TO LWB a DO DIGIT carry := 0;
- calculate each digit (whilst removing the carry) #
FOR place b FROM UPB b BY next digit TO LWB b DO
# n.b. result may be 2 digits #
INT result := ab[place a + place b] + a[place a]*b[place b] + carry;
carry := result % digit base; # avoid MOD as it doesn't under handle -ve numbers #
ab[place a + place b] := result - carry * digit base
OD;
ab[place a + LWB b + next digit] +:= carry
OD; NORMALISE ab
);</lang>
<lang algol68># The following standard operators could (potentially) also be defined # OP - = (INTEGER a)INTEGER: raise integer not implemented error("monadic minus"),
ABS = (INTEGER a)INTEGER: raise integer not implemented error("ABS"),
ODD = (INTEGER a)INTEGER: raise integer not implemented error("ODD"),
BIN = (INTEGER a)INTEGER: raise integer not implemented error("BIN");
OP + = (INTEGER a, b)INTEGER: raise integer not implemented error("addition"),
- = (INTEGER a, b)INTEGER: raise integer not implemented error("subtraction"),
/ = (INTEGER a, b)REAL: ( VOID(raise integer not implemented error("floating point division")); SKIP),
% = (INTEGER a, b)INTEGER: raise integer not implemented error("fixed point division"),
%* = (INTEGER a, b)INTEGER: raise integer not implemented error("modulo division"),
** = (INTEGER a, b)INTEGER: raise integer not implemented error("to the power of");
LONG INT default integer width := long long int width - 2;
INT fix w = -1177584, fix h = -3915074; # floating point error, probably GMP/hardware specific #
INTEGER golden ratio w := INTEGERINIT ENTIER ((long long sqrt(5)-1) / 2 * LENG LENG 10 ** default integer width + fix w),
golden ratio h := INTEGERINIT ENTIER ((long long sqrt(5)+1) / 2 * LENG LENG 10 ** default integer width + fix h);
test: (
print(( "The approximate golden ratios, width: ", REPR golden ratio w, new line, " length: ", REPR golden ratio h, new line, " product is exactly: ", REPR (golden ratio w * golden ratio h), new line));
INTEGER two to the power of 64 = INTEGERINIT(LONG 2 ** 64);
INTEGER neg two to the power of 64 = INTEGERINIT(-(LONG 2 ** 64));
print(("2 ** 64 * -(2 ** 64) = ", REPR (two to the power of 64 * neg two to the power of 64), new line))
)</lang>
Output:
The approximate golden ratios, width: +618,033,988,749,894,848,204,586,834,365,638,117,720,309,179,805,762,862,135,448,622,705,261
length: +1,618,033,988,749,894,848,204,586,834,365,638,117,720,309,179,805,762,862,135,448,622,705,261
product is exactly: +1,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,001,201,173,450,350,400,438,606,015,942,314,498,798,603,569,682,901,026,716,145,698,077,078,121
2 ** 64 * -(2 ** 64) = -340,282,366,920,938,463,463,374,607,431,768,211,456
Other libraries or implementation specific extensions
As of February 2009 no open source libraries to do this task have been located.
AutoHotkey
ahk discussion <lang autohotkey>MsgBox % x := mul(256,256) MsgBox % x := mul(x,x) MsgBox % x := mul(x,x) ; 18446744073709551616 MsgBox % x := mul(x,x) ; 340282366920938463463374607431768211456
mul(b,c) { ; <- b*c
VarSetCapacity(a, n:=StrLen(b)+StrLen(c), 48), NumPut(0,a,n,"char")
Loop % StrLen(c) {
i := StrLen(c)+1-A_Index, cy := 0
Loop % StrLen(b) {
j := StrLen(b)+1-A_Index,
t := SubStr(a,i+j,1) + SubStr(b,j,1) * SubStr(c,i,1) + cy
cy := t // 10
NumPut(mod(t,10)+48,a,i+j-1,"char")
}
NumPut(cy+48,a,i+j-2,"char")
}
Return cy ? a : SubStr(a,2)
}</lang>
AWK
<lang awk>BEGIN {
DEBUG = 0
n = 2^64
nn = sprintf("%.0f", n)
printf "2^64 * 2^64 = %.0f\n", multiply(nn, nn)
printf "2^64 * 2^64 = %.0f\n", n*n
exit
}
function multiply(x, y, len_x,len_y,ax,ay,j,m,c,i,k,d,v,res,mul,result) {
len_x = split_reverse(x, ax)
len_y = split_reverse(y, ay)
print_array(ax)
print_array(ay)
for (j=1; j<=len_y; j++) {
m = ay[j]
c = 0
i = j - 1
for (k=1; k<=len_x; k++) {
d = ax[k]
i++
v = res[i]
if (v == "") {
append_array(res, 0)
v = 0
}
mul = v + c + d*m
c = int(mul / 10)
v = mul % 10
res[i] = v
}
append_array(res, c)
}
print_array(res)
result = reverse_join(res)
sub(/^0+/, "", result)
return result
}
function split_reverse(x, a, a_x) {
split(x, a_x, "") return reverse_array(a_x, a)
}
function reverse_array(a,b, len,i) {
len = length_array(a)
for (i in a) {
b[1+len-i] = a[i]
}
return len
}
function length_array(a, len,i) {
len = 0 for (i in a) len++ return len
}
function append_array(a, value, len) {
len = length_array(a) a[++len] = value
}
function reverse_join(a, len,str,i) {
len = length_array(a)
str = ""
for (i=len; i>=1; i--) {
str = str a[i]
}
return str
}
function print_array(a, len,i) {
if (DEBUG) {
len = length_array(a)
print "length=" len
for (i=1; i<=len; i++) {
printf("%s ", i%10)
}
print ""
for (i=1; i<=len; i++) {
#print i " " a[i]
printf("%s ", a[i])
}
print ""
print "===="
}
}</lang> outputs:
2^64 * 2^64 = 340282366920938463463374607431768211456 2^64 * 2^64 = 340282366920938463463374607431768211456
BASIC
Version 1
<lang qbasic>'PROGRAM : BIG MULTIPLICATION VER #1 'LRCVS 01.01.2010 'THIS PROGRAM SIMPLY MAKES A MULTIPLICATION 'WITH ALL THE PARTIAL PRODUCTS. '............................................................
DECLARE SUB A.INICIO (A$, B$) DECLARE SUB B.STORE (CAD$, N$) DECLARE SUB C.PIZARRA () DECLARE SUB D.ENCABEZADOS (A$, B$) DECLARE SUB E.MULTIPLICACION (A$, B$) DECLARE SUB G.SUMA () DECLARE FUNCTION F.INVCAD$ (CAD$)
RANDOMIZE TIMER CALL A.INICIO(A$, B$) CALL B.STORE(A$, "A") CALL B.STORE(B$, "B") CALL C.PIZARRA CALL D.ENCABEZADOS(A$, B$) CALL E.MULTIPLICACION(A$, B$) CALL G.SUMA
SUB A.INICIO (A$, B$)
CLS
'Note: Number of digits > 1000
INPUT "NUMBER OF DIGITS "; S
CLS
A$ = ""
B$ = ""
FOR N = 1 TO S
A$ = A$ + LTRIM$(STR$(INT(RND * 9)))
NEXT N
FOR N = 1 TO S
B$ = B$ + LTRIM$(STR$(INT(RND * 9)))
NEXT N
END SUB
SUB B.STORE (CAD$, N$)
OPEN "O", #1, N$
FOR M = LEN(CAD$) TO 1 STEP -1
WRITE #1, MID$(CAD$, M, 1)
NEXT M
CLOSE (1)
END SUB
SUB C.PIZARRA
OPEN "A", #3, "R" WRITE #3, "" CLOSE (3) KILL "R"
END SUB
SUB D.ENCABEZADOS (A$, B$)
LT = LEN(A$) + LEN(B$) + 1 L$ = STRING$(LT, " ") OPEN "A", #3, "R" MID$(L$, LT - LEN(A$) + 1) = A$ WRITE #3, L$ CLOSE (3) L$ = STRING$(LT, " ") OPEN "A", #3, "R" MID$(L$, LT - LEN(B$) - 1) = "X " + B$ WRITE #3, L$ CLOSE (3)
END SUB
SUB E.MULTIPLICACION (A$, B$)
LT = LEN(A$) + LEN(B$) + 1
L$ = STRING$(LT, " ")
C$ = ""
D$ = ""
E$ = ""
CT1 = 1
ACUM = 0
OPEN "I", #2, "B"
WHILE EOF(2) <> -1
INPUT #2, B$
OPEN "I", #1, "A"
WHILE EOF(1) <> -1
INPUT #1, A$
RP = (VAL(A$) * VAL(B$)) + ACUM
C$ = LTRIM$(STR$(RP))
IF EOF(1) <> -1 THEN D$ = D$ + RIGHT$(C$, 1)
IF EOF(1) = -1 THEN D$ = D$ + F.INVCAD$(C$)
E$ = LEFT$(C$, LEN(C$) - 1)
ACUM = VAL(E$)
WEND
CLOSE (1)
MID$(L$, LT - CT1 - LEN(D$) + 2) = F.INVCAD$(D$)
OPEN "A", #3, "R"
WRITE #3, L$
CLOSE (3)
L$ = STRING$(LT, " ")
ACUM = 0
C$ = ""
D$ = ""
E$ = ""
CT1 = CT1 + 1
WEND
CLOSE (2)
END SUB
FUNCTION F.INVCAD$ (CAD$)
LCAD = LEN(CAD$)
CADTEM$ = ""
FOR CAD = LCAD TO 1 STEP -1
CADTEM$ = CADTEM$ + MID$(CAD$, CAD, 1)
NEXT CAD
F.INVCAD$ = CADTEM$
END FUNCTION
SUB G.SUMA
CF = 0
OPEN "I", #3, "R"
WHILE EOF(3) <> -1
INPUT #3, R$
CF = CF + 1
AN = LEN(R$)
WEND
CF = CF - 2
CLOSE (3)
W$ = ""
ST = 0
ACUS = 0
FOR P = 1 TO AN
K = 0
OPEN "I", #3, "R"
WHILE EOF(3) <> -1
INPUT #3, R$
K = K + 1
IF K > 2 THEN ST = ST + VAL(MID$(R$, AN - P + 1, 1))
IF K > 2 THEN M$ = LTRIM$(STR$(ST + ACUS))
WEND
'COLOR 10: LOCATE CF + 3, AN - P + 1: PRINT RIGHT$(M$, 1); : COLOR 7
W$ = W$ + RIGHT$(M$, 1)
ACUS = VAL(LEFT$(M$, LEN(M$) - 1))
CLOSE (3)
ST = 0
NEXT P
OPEN "A", #3, "R" WRITE #3, " " + RIGHT$(F.INVCAD(W$), AN - 1) CLOSE (3) CLS PRINT "THE SOLUTION IN THE FILE: R"
END SUB</lang>
Version 2
<lang qbasic>'PROGRAM: BIG MULTIPLICATION VER # 2 'LRCVS 01/01/2010 'THIS PROGRAM SIMPLY MAKES A BIG MULTIPLICATION 'WITHOUT THE PARTIAL PRODUCTS. 'HERE SEE ONLY THE SOLUTION. '............................................................... CLS PRINT "WAIT"
NA = 2000 'NUMBER OF ELEMENTS OF THE MULTIPLY. NB = 2000 'NUMBER OF ELEMENTS OF THE MULTIPLIER. 'Solution = 4000 Exacts digits
'...................................................... OPEN "X" + ".MLT" FOR BINARY AS #1 CLOSE (1) KILL "*.MLT" '..................................................... 'CREATING THE MULTIPLY >>> A 'CREATING THE MULTIPLIER >>> B FOR N = 1 TO 2 IF N = 1 THEN F$ = "A" + ".MLT": NN = NA IF N = 2 THEN F$ = "B" + ".MLT": NN = NB
OPEN F$ FOR BINARY AS #1
FOR N2 = 1 TO NN
RANDOMIZE TIMER
X$ = LTRIM$(STR$(INT(RND * 10)))
SEEK #1, N2: PUT #1, N2, X$
NEXT N2
SEEK #1, N2
CLOSE (1)
NEXT N '..................................................... OPEN "A" + ".MLT" FOR BINARY AS #1 FOR K = 0 TO 9 NUM$ = "": Z$ = "": ACU = 0: GG = NA C$ = LTRIM$(STR$(K))
OPEN C$ + ".MLT" FOR BINARY AS #2
'OPEN "A" + ".MLT" FOR BINARY AS #1
FOR N = 1 TO NA
SEEK #1, GG: GET #1, GG, X$
NUM$ = X$
Z$ = LTRIM$(STR$(ACU + (VAL(X$) * VAL(C$))))
L = LEN(Z$)
ACU = 0
IF L = 1 THEN NUM$ = Z$: PUT #2, N, NUM$
IF L > 1 THEN ACU = VAL(LEFT$(Z$, LEN(Z$) - 1)): NUM$ = RIGHT$(Z$, 1): PUT #2, N, NUM$
SEEK #2, N: PUT #2, N, NUM$
GG = GG - 1
NEXT N
IF L > 1 THEN ACU = VAL(LEFT$(Z$, LEN(Z$) - 1)): NUM$ = LTRIM$(STR$(ACU)): XX$ = XX$ + NUM$: PUT #2, N, NUM$
'CLOSE (1)
CLOSE (2)
NEXT K CLOSE (1) '...................................................... ACU = 0 LT5 = 1 LT6 = LT5 OPEN "B" + ".MLT" FOR BINARY AS #1
OPEN "D" + ".MLT" FOR BINARY AS #3
FOR JB = NB TO 1 STEP -1
SEEK #1, JB
GET #1, JB, X$
OPEN X$ + ".MLT" FOR BINARY AS #2: LF = LOF(2): CLOSE (2)
OPEN X$ + ".MLT" FOR BINARY AS #2
FOR KB = 1 TO LF
SEEK #2, KB
GET #2, , NUM$
SEEK #3, LT5
GET #3, LT5, PR$
T$ = ""
T$ = LTRIM$(STR$(ACU + VAL(NUM$) + VAL(PR$)))
PR$ = RIGHT$(T$, 1)
ACU = 0
IF LEN(T$) > 1 THEN ACU = VAL(LEFT$(T$, LEN(T$) - 1))
SEEK #3, LT5: PUT #3, LT5, PR$
LT5 = LT5 + 1
NEXT KB
IF ACU <> 0 THEN PR$ = LTRIM$(STR$(ACU)): PUT #3, LT5, PR$
CLOSE (2)
LT6 = LT6 + 1
LT5 = LT6
ACU = 0
NEXT JB
CLOSE (3)
CLOSE (1) OPEN "D" + ".MLT" FOR BINARY AS #3: LD = LOF(3): CLOSE (3) ER = 1 OPEN "D" + ".MLT" FOR BINARY AS #3
OPEN "R" + ".MLT" FOR BINARY AS #4
FOR N = LD TO 1 STEP -1
SEEK #3, N: GET #3, N, PR$
SEEK #4, ER: PUT #4, ER, PR$
ER = ER + 1
NEXT N
CLOSE (4)
CLOSE (3) KILL "D.MLT" FOR N = 0 TO 9
C$ = LTRIM$(STR$(N)) KILL C$ + ".MLT"
NEXT N PRINT "END" PRINT "THE SOLUTION IN THE FILE: R.MLT"</lang>
C
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <string.h>
- define X(V,L,I) ( ((I)<(L)) ? (V[(L)-(I)-1]-'0') : 0)
- define MIN(A,B) ( ((A)<(B)) ? (A) : (B) )
char *longmult(const char *a, const char *b) {
int n, m, T; char *r = NULL; int i, j, k; int *C, *R; n = strlen(a); m = strlen(b); T = n+m+2;
r = malloc(T+1); R = malloc(T*sizeof(int)); C = malloc((n+1)*(m+1)*sizeof(int));
memset(r, '0', T); memset(C, 0, (n+1)*(m+1)*sizeof(int)); r[T] = 0;
for(i=0; i<(m+1); i++)
{
C[i*(n+1)] = (X(b,m,i) * X(a,n,0));
for(j=1; j<(n+1); j++)
{
C[i*(n+1)+j] = (X(b,m,i) * X(a,n,j) + C[i*(n+1)+j-1] / 10);
}
}
for(k=0; k < T; k++)
{
R[k] = 0;
for(j=0; j < MIN(k,m+1) ; j++)
{
i = k-j-1;
if ( (i>n) || (i<0) ) { continue; }
R[k] += (C[j*(n+1)+i]%10);
}
R[k] += ((k-1)<0 ) ? 0 : (R[k-1]/10);
}
for(k=T; k>0; k--) r[k-1] = R[T-k+1]%10 + '0';
free(C); free(R);
return r;
}</lang>
<lang c>/* using */ const char *n1 = "18446744073709551616"; int main() {
char *res; int lz;
/* printf("%s * %s = 340282366920938463463374607431768211456\n", n1, n1); */
res = longmult(n1,n1);
for(lz=0; (lz < strlen(res)) && (res[lz]=='0') ;lz++) ;
printf("%s * %s = %s\n", n1, n1, res+lz);
free(res);
return 0;
}</lang>
The code does not handle negative intergers, nor there are proper error checks. It is more or less the implementation of the way a human being multiplies two integers. The numbers are stored as strings.
Using GMP (GNU Multi Precision library)
<lang c>#include <stdio.h>
- include <gmp.h>
int main() {
mpz_t z1, z2, zr;
mpz_init(z2); mpz_init(zr);
mpz_init_set_str(z1, "18446744073709551616", 10);
mpz_set(z2, z1);
mpz_mul(zr, z1, z2);
mpz_out_str(stdout, 10, zr);
printf("\n");
return 0;
}</lang>
C++
using the Class Library for Numbers , CLN, and g++ <lang cpp>#include <cln/integer.h>//for mathematical operations on arbitrarily long integers
- include <cln/integer_io.h>//for input/output of long integers
- include <iostream>
int main( ) {
cln::cl_I base = 2 , exponent = 64 ;//cln is a namespace cln::cl_I factor = cln::expt_pos( base , exponent ) ; cln::cl_I product = factor * factor ; std::cout << "The result of 2^64 * 2^64 is " << product << " !\n" ; return 0 ;
}</lang>
Output: The result of 2^64 * 2^64 is 340282366920938463463374607431768211456 !
C#
System.Numerics.BigInteger was added with C# 4.
<lang csharp>using System; using System.Numerics;
class Program {
static void Main() {
BigInteger pow2_64 = BigInteger.Pow(2, 64);
BigInteger result = BigInteger.Multiply(pow2_64, pow2_64);
Console.WriteLine(result);
}
}</lang>
Output:
340282366920938463463374607431768211456
Common Lisp
<lang lisp>(defun number->digits (number)
(do ((digits '())) ((zerop number) digits)
(multiple-value-bind (quotient remainder) (floor number 10)
(setf number quotient)
(push remainder digits))))
(defun digits->number (digits)
(reduce #'(lambda (n d) (+ (* 10 n) d)) digits :initial-value 0))
(defun long-multiply (a b)
(labels ((first-digit (list)
"0 if list is empty, else first element of list."
(if (endp list) 0
(first list)))
(long-add (digitses &optional (carry 0) (sum '()))
"Do long addition on the list of lists of digits. Each
list of digits in digitses should begin with the least
significant digit. This is the opposite of the digit
list returned by number->digits which places the most
significant digit first. The digits returned by
long-add do have the most significant bit first."
(if (every 'endp digitses)
(nconc (digits carry) sum)
(let ((column-sum (reduce '+ (mapcar #'first-digit digitses)
:initial-value carry)))
(multiple-value-bind (carry column-digit)
(floor column-sum 10)
(long-add (mapcar 'rest digitses)
carry (list* column-digit sum)))))))
;; get the digits of a and b (least significant bit first), and
;; compute the zero padded rows. Then, add these rows (using
;; long-add) and convert the digits back to a number.
(do ((a (nreverse (digits a)))
(b (nreverse (digits b)))
(prefix '() (list* 0 prefix))
(rows '()))
((endp b) (digits->number (long-add rows)))
(let* ((bi (pop b))
(row (mapcar #'(lambda (ai) (* ai bi)) a)))
(push (append prefix row) rows)))))</lang>
> (long-multiply (expt 2 64) (expt 2 64)) 340282366920938463463374607431768211456
D
Using the standard library: <lang d>import std.stdio, std.bigint;
void main() {
writeln(BigInt(2) ^^ 64 * BigInt(2) ^^ 64);
}</lang> Output:
340282366920938463463374607431768211456
Long multiplication, same output:
<lang d>import std.stdio, std.algorithm, std.conv, std.range;
auto longMult(string x, string y) {
auto digits1 = array(map!q{a - '0'}(retro(x)));
auto digits2 = array(map!q{a - '0'}(retro(y)));
int[] res;
foreach (i, d1; digits1)
foreach (j, d2; digits2) {
int idx = i + j;
if (res.length <= idx)
res.length += 1;
res[idx] = d1 * d2 + res[idx];
if (res[idx] > 9) {
if (res.length <= idx + 1)
res.length += 1;
res[idx + 1] = res[idx] / 10 + res[idx + 1];
res[idx] -= res[idx] / 10 * 10;
}
}
return map!q{ cast(char)(a + '0') }(retro(res));
}
void main() {
string two64 = "18446744073709551616"; writeln(longMult(two64, two64));
}</lang>
Euphoria
<lang euphoria>constant base = 1000000000
function atom_to_long(atom a)
sequence s
s = {}
while a>0 do
s = append(s,remainder(a,base))
a = floor(a/base)
end while
return s
end function
function long_mult(object a, object b)
sequence c
if atom(a) then
a = atom_to_long(a)
end if
if atom(b) then
b = atom_to_long(b)
end if
c = repeat(0,length(a)+length(b))
for i = 1 to length(a) do
c[i .. i+length(b)-1] += a[i]*b
end for
for i = 1 to length(c) do
if c[i] > base then
c[i+1] += floor(c[i]/base) -- carry
c[i] = remainder(c[i],base)
end if
end for
if c[$] = 0 then
c = c[1..$-1]
end if
return c
end function
function long_to_str(sequence a)
sequence s
s = sprintf("%d",a[$])
for i = length(a)-1 to 1 by -1 do
s &= sprintf("%09d",a[i])
end for
return s
end function
sequence a, b, c
a = atom_to_long(power(2,32)) printf(1,"a is %s\n",{long_to_str(a)})
b = long_mult(a,a) printf(1,"a*a is %s\n",{long_to_str(b)})
c = long_mult(b,b) printf(1,"a*a*a*a is %s\n",{long_to_str(c)})</lang>
Output:
a is 4294967296 a*a is 18446744073709551616 a*a*a*a is 340282366920938463488374607424768211456
F#
<lang F#>> let X = 2I ** 64 * 2I ** 64 ;;
val X : System.Numerics.BigInteger = 340282366920938463463374607431768211456 </lang>
Fortran
<lang fortran>module LongMoltiplication
implicit none
type longnum
integer, dimension(:), pointer :: num
end type longnum
interface operator (*)
module procedure longmolt_ll
end interface
contains
subroutine longmolt_s2l(istring, num) character(len=*), intent(in) :: istring type(longnum), intent(out) :: num integer :: i, l
l = len(istring)
allocate(num%num(l))
forall(i=1:l) num%num(l-i+1) = iachar(istring(i:i)) - 48
end subroutine longmolt_s2l
! this one performs the moltiplication function longmolt_ll(a, b) result(c) type(longnum) :: c type(longnum), intent(in) :: a, b integer, dimension(:,:), allocatable :: t integer :: ntlen, i, j
ntlen = size(a%num) + size(b%num) + 1 allocate(c%num(ntlen)) c%num = 0
allocate(t(size(b%num), ntlen)) t = 0 forall(i=1:size(b%num), j=1:size(a%num)) t(i, j+i-1) = b%num(i) * a%num(j)
do j=2, ntlen
forall(i=1:size(b%num)) t(i, j) = t(i, j) + t(i, j-1)/10
end do
forall(j=1:ntlen) c%num(j) = sum(mod(t(:,j), 10))
do j=2, ntlen
c%num(j) = c%num(j) + c%num(j-1)/10
end do
c%num = mod(c%num, 10) deallocate(t) end function longmolt_ll
subroutine longmolt_print(num) type(longnum), intent(in) :: num
integer :: i, j
do j=size(num%num), 2, -1
if ( num%num(j) /= 0 ) exit
end do
do i=j, 1, -1
write(*,"(I1)", advance="no") num%num(i)
end do
end subroutine longmolt_print
end module LongMoltiplication</lang>
<lang fortran>program Test
use LongMoltiplication
type(longnum) :: a, b, r
call longmolt_s2l("18446744073709551616", a)
call longmolt_s2l("18446744073709551616", b)
r = a * b call longmolt_print(r) write(*,*)
end program Test</lang>
Haskell
<lang haskell>digits :: Integer -> [Integer] digits = map (fromIntegral.digitToInt) . show
lZZ = inits $ repeat 0
table f = map . flip (map . f)
polymul = ((map sum . transpose . zipWith (++) lZZ) .) . table (*)
longmult = (foldl1 ((+) . (10 *)) .) . (. digits) . polymul . digits</lang> Output: <lang haskell>*Main> (2^64) `longmult` (2^64) 340282366920938463463374607431768211456</lang>
Icon and Unicon
Icon and Unicon support large integers. <lang Icon>procedure main() write(2^64*2^64) end</lang>
J
Solution: <lang j> digits =: ,.&.":
polymult =: +//.@(*/) buildDecimal=: 10x&#.
longmult=: buildDecimal@polymult&digits</lang>
Example: <lang j> longmult~ 2x^64 340282366920938463463374607431768211456</lang>
Alternatives:
longmult could have been defined concisely:
<lang j>longmult=: 10x&#.@(+//.@(*/)&(,.&.":))</lang>
Or, of course, the task may be accomplished without the verb definitions:
<lang j> 10x&#.@(+//.@(*/)&(,.&.":))~2x^64
340282366920938463463374607431768211456</lang>
Or using the code (+ 10x&*)/@|. instead of #.:
<lang j> (+ 10x&*)/@|.@(+//.@(*/)&(,.&.":))~2x^64
340282366920938463463374607431768211456</lang>
Or you could use the built-in language support for arbitrary precision multiplication:
<lang j> (2x^64)*(2x^64)
340282366920938463463374607431768211456</lang>
Explaining the component verbs:
digitstranslates a number to a corresponding list of digits;
<lang j> ,.&.": 123 1 2 3</lang>
polymult(multiplies polynomials): ref. [1]
<lang j> 1 2 3 (+//.@(*/)) 1 2 3 1 4 10 12 9</lang>
buildDecimal(translates a list of decimal digits to the corresponding extended precision number):
<lang j> (+ 10x&*)/|. 1 4 10 12 9 15129</lang>
Java
This is a straight-forward implementation of Long multiplication. It works with numbers of any length since it uses BigInteger. <lang java>import java.math.BigInteger;
public class LongMult {
public static void main(String[] args) { BigInteger TwoPow64 = new BigInteger("18446744073709551616"); System.out.println(mult(TwoPow64, TwoPow64)); }
public static BigInteger mult(BigInteger a, BigInteger b){ return a.multiply(b); } }</lang> Output:
340282366920938463463374607431768211456
JavaScript
<lang javascript>function mult(num1,num2){ var a1 = num1.split("").reverse(); var a2 = num2.split("").reverse(); var aResult = new Array;
for ( iterNum1 = 0; iterNum1 < a1.length; iterNum1++ ) { for ( iterNum2 = 0; iterNum2 < a2.length; iterNum2++ ) { idxIter = iterNum1 + iterNum2; // Get the current array position. aResult[idxIter] = a1[iterNum1] * a2[iterNum2] + ( idxIter >= aResult.length ? 0 : aResult[idxIter] );
if ( aResult[idxIter] > 9 ) { // Carrying aResult[idxIter + 1] = Math.floor( aResult[idxIter] / 10 ) + ( idxIter + 1 >= aResult.length ? 0 : aResult[idxIter + 1] ); aResult[idxIter] -= Math.floor( aResult[idxIter] / 10 ) * 10; } } } return aResult.reverse().join(""); }</lang>
Liberty BASIC
<lang lb>print 2^64 *2^64</lang>
Mathematica
We define the long multiplication function: <lang Mathematica> LongMultiplication[a_,b_]:=Module[{d1,d2},
d1=IntegerDigits[a]//Reverse; d2=IntegerDigits[b]//Reverse; Sum[d1id2j*10^(i+j-2),{i,1,Length[d1]},{j,1,Length[d2]}] ]</lang>
Example: <lang Mathematica> n1 = 2^64;
n2 = 2^64; LongMultiplication[n1, n2]</lang>
gives back: <lang Mathematica> 340282366920938463463374607431768211456</lang>
To check the speed difference between built-in multiplication (which is already arbitrary precision) we multiply two big numbers (2^8000 has 2409 digits!) and divide their timings: <lang Mathematica> n1=2^8000;
n2=2^8000; Timing[LongMultiplication[n1,n2]]1 Timing[n1 n2]1 Floor[%%/%]</lang>
gives back: <lang Mathematica> 72.9686
7.*10^-6 10424088</lang>
So our custom function takes about 73 second, the built-in function a couple of millionths of a second, so the long multiplication is about 10.5 million times slower! Mathematica uses Karatsuba multiplication for large integers, which is several magnitudes faster for really big numbers. Making it able to multiply in about a second; the final result has 9542426 digits; result omitted for obvious reasons.
Perl
<lang perl>#!/usr/bin/perl -w use strict;
- This should probably be done in a loop rather than be recursive.
sub add_with_carry {
my $resultref = shift; my $addend = shift; my $addendpos = shift;
push @$resultref, (0) while (scalar @$resultref < $addendpos + 1); my $addend_result = $addend + $resultref->[$addendpos]; my @addend_digits = reverse split //, $addend_result; $resultref->[$addendpos] = shift @addend_digits;
my $carry_digit = shift @addend_digits; &add_with_carry($resultref, $carry_digit, $addendpos + 1) if( defined $carry_digit )
}
sub longhand_multiplication {
my @multiplicand = reverse split //, shift;
my @multiplier = reverse split //, shift;
my @result = ();
my $multiplicand_offset = 0;
foreach my $multiplicand_digit (@multiplicand)
{
my $multiplier_offset = $multiplicand_offset;
foreach my $multiplier_digit (@multiplier)
{
my $multiplication_result = $multiplicand_digit * $multiplier_digit;
my @result_digit_addend_list = reverse split //, $multiplication_result;
my $addend_offset = $multiplier_offset;
foreach my $result_digit_addend (@result_digit_addend_list)
{
&add_with_carry(\@result, $result_digit_addend, $addend_offset++)
}
++$multiplier_offset; }
++$multiplicand_offset; }
@result = reverse @result;
return join , @result;
}
my $sixtyfour = "18446744073709551616";
my $onetwentyeight = &longhand_multiplication($sixtyfour, $sixtyfour); print "$onetwentyeight\n";</lang>
Perl 6
For efficiency (and novelty), this program explicitly implements long multiplication, but in base 10000. That base was chosen because multiplying two 5-digit numbers can overflow a 32-bit integer, but two 4-digit numbers cannot. <lang perl6>sub num_to_groups ( $num ) { $num.flip.comb(/.**1..4/)».flip }; sub groups_to_num ( @g ) { [~] @g.pop, @g.reverse».fmt('%04d') };
sub long_multiply ( Str $x, Str $y ) {
my @group_sums;
for num_to_groups($x).pairs X num_to_groups($y).pairs -> $xp, $yp {
@group_sums[ $xp.key + $yp.key ] += $xp.value * $yp.value;
}
for @group_sums.keys -> $k {
next if @group_sums[$k] < 10000;
@group_sums[$k+1] += @group_sums[$k].Int div 10000;
@group_sums[$k] %= 10000;
}
return groups_to_num @group_sums;
}
long_multiply( '18446744073709551616', '18446744073709551616' ).say;</lang>
In any case, integers are specced to be arbitrarily large, which at least one implementation supports already:
pugs> 18446744073709551616 * 18446744073709551616 340282366920938463463374607431768211456
PHP
<lang PHP><?php
$factor = bcpow(2, 64); $product = bcmul($factor, $factor); echo "2^64 * 2^64 is " . $product; ?></lang> Output:
2^64 * 2^64 is 340282366920938463463374607431768211456
PicoLisp
<lang PicoLisp>: (* (** 2 64) (** 2 64)) -> 340282366920938463463374607431768211456</lang>
PL/I
<lang PL/I>/* Multiply a by b, giving c. */ multiply: procedure (a, b, c);
declare (a, b, c) (*) fixed decimal (1); declare (d, e, f) (hbound(a,1)) fixed decimal (1); declare pr (-hbound(a,1) : hbound(a,1)) fixed decimal (1); declare p fixed decimal (2), (carry, s) fixed decimal (1); declare neg bit (1) aligned; declare (i, j, n, offset) fixed binary (31);
n = hbound(a,1);
d = a;
e = b;
s = a(1) + b(1);
neg = (s = 9);
if a(1) = 9 then call complement (d);
if b(1) = 9 then call complement (e);
pr = 0;
offset = 0; carry = 0;
do i = n to 1 by -1;
do j = n to 1 by -1;
p = d(i) * e(j) + pr(j-offset) + carry;
if p > 9 then do; carry = p/10; p = mod(p, 10); end; else carry = 0;
pr(j-offset) = p;
end;
offset = offset + 1;
end;
do i = hbound(a,1) to 1 by -1;
c(i) = pr(i);
end;
do i = -hbound(a,1) to 1;
if pr(i) ^= 0 then signal fixedoverflow;
end;
if neg then call complement (c);
end multiply;
complement: procedure (a);
declare a(*) fixed decimal (1); declare i fixed binary (31), carry fixed decimal (1); declare s fixed decimal (2);
carry = 1;
do i = hbound(a,1) to 1 by -1;
s = 9 - a(i) + carry;
if s > 9 then do; s = s - 10; carry = 1; end; else carry = 0;
a(i) = s;
end;
end complement;</lang> Calling sequence: <lang PL/I> a = 0; b = 0; c = 0;
a(60) = 1;
do i = 1 to 64; /* Generate 2**64 */
call add (a, a, b);
put skip;
call output (b);
a = b;
end;
call multiply (a, b, c);
put skip;
call output (c);</lang>
Final output:
18446744073709551616 340282366920938463463374607431768211456
PowerShell
<lang PowerShell># LongAddition only supports Unsigned Integers represented as Strings/Character Arrays Function LongAddition ( [Char[]] $lhs, [Char[]] $rhs ) { $lhsl = $lhs.length $rhsl = $rhs.length if(($lhsl -gt 0) -and ($rhsl -gt 0)) { $maxplace = [Math]::Max($rhsl,$lhsl)+1 1..$maxplace | ForEach-Object { $carry = 0 $result = "" } { $add1 = 0 $add2 = 0 if( $_ -le $lhsl ) { $add1 = [int]$lhs[ -$_ ] - 48 } if( $_ -le $rhsl ) { $add2 = [int]$rhs[ -$_ ] - 48 } $iresult = $add1 + $add2 + $carry if( ( $_ -lt $maxplace ) -or ( $iresult -gt 0 ) ) { $result = "{0}{1}" -f ( $iresult % 10 ),$result } $carry = [Math]::Floor( $iresult / 10 ) } { $result } } elseif($lhsl -gt 0) { [String]::Join( , $lhs ) } elseif($rhsl -gt 0) { [String]::Join( , $rhs ) } else { "0" } }
- LongMultiplication only supports Unsigned Integers represented as Strings/Character Arrays
Function LongMultiplication ( [Char[]] $lhs, [Char[]] $rhs ) { $lhsl = $lhs.length $rhsl = $rhs.length if(($lhsl -gt 0) -and ($rhsl -gt 0)) { 1..$lhsl | ForEach-Object { $carry0 = "" $result0 = "" } { $i = -$_ $add1 = ( 1..$rhsl | ForEach-Object { $carry1 = 0 $result1 = "" } { $j = -$_ $mult1 = [int]$lhs[ $i ] - 48 $mult2 = [int]$rhs[ $j ] - 48 $iresult1 = $mult1 * $mult2 + $carry1 $result1 = "{0}{1}" -f ( $iresult1 % 10 ), $result1 $carry1 = [Math]::Floor( $iresult1 / 10 ) } { if( $carry1 -gt 0 ) { $result1 = "{0}{1}" -f $carry1, $result1 } $result1 } ) $iresult0 = ( LongAddition $add1 $carry0 ) $iresultl = $iresult0.length $result0 = "{0}{1}" -f $iresult0[-1],$result0 if( $iresultl -gt 1 ) { $carry0 = [String]::Join( , $iresult0[ -$iresultl..-2 ] ) } else { $carry0 = "" } } { if( $carry0 -ne "" ) { $result0 = "{0}{1}" -f $carry0, $result0 } $result0 } } else { "0" } }
LongMultiplication "18446744073709551616" "18446744073709551616"</lang>
Prolog
<lang Prolog> ?- X is 2**64 * 2**64. X = 340282366920938463463374607431768211456.</lang>
PureBasic
Using Decimal.pbi by Stargåte allows for calculation with long numbers, this is useful since version 4.41 of PureBasic mostly only supporter data types native to x86/x64/PPC etc processors. <lang PureBasic>XIncludeFile "decimal.pbi"
Define.Decimal *a, *b
- a=PowerDecimal(IntegerToDecimal(2),IntegerToDecimal(64))
- b=TimesDecimal(*a,*a,#NoDecimal)
Print("2^64*2^64 = "+DecimalToString(*b))</lang>
Outputs
2^64*2^64 = 340282366920938463463374607431768211456
Python
(Note that Python comes with arbitrary length integers).
<lang python>#!/usr/bin/env python print 2**64*2**64</lang>
<lang python>#!/usr/bin/env python
def add_with_carry(result, addend, addendpos):
while True:
while len(result) < addendpos + 1:
result.append(0)
addend_result = str(int(addend) + int(result[addendpos]))
addend_digits = list(addend_result)
result[addendpos] = addend_digits.pop()
if not addend_digits:
break
addend = addend_digits.pop()
addendpos += 1
def longhand_multiplication(multiplicand, multiplier):
result = []
for multiplicand_offset, multiplicand_digit in enumerate(reversed(multiplicand)):
for multiplier_offset, multiplier_digit in enumerate(reversed(multiplier), start=multiplicand_offset):
multiplication_result = str(int(multiplicand_digit) * int(multiplier_digit))
for addend_offset, result_digit_addend in enumerate(reversed(multiplication_result), start=multiplier_offset):
add_with_carry(result, result_digit_addend, addend_offset)
result.reverse()
return .join(result)
if __name__ == "__main__":
sixtyfour = "18446744073709551616"
onetwentyeight = longhand_multiplication(sixtyfour, sixtyfour) print(onetwentyeight)</lang>
Shorter version:
<lang python>#!/usr/bin/env python
def digits(x):
return [int(c) for c in str(x)]
def mult_table(xs, ys):
return [[x * y for x in xs] for y in ys]
def polymul(xs, ys):
return map(lambda *vs: sum(filter(None, vs)),
*[[0] * i + zs for i, zs in enumerate(mult_table(xs, ys))])
def longmult(x, y):
result = 0
for v in polymul(digits(x), digits(y)):
result = result * 10 + v
return result
if __name__ == "__main__":
print longmult(2**64, 2**64)</lang>
R
Using GMP
<lang R>library(gmp) a <- as.bigz("18446744073709551616") mul.bigz(a,a)</lang>
"340282366920938463463374607431768211456"
A native implementation
This code is more verbose than necessary, for ease of understanding. <lang R>longmult <- function(xstr, ystr) {
#get the number described in each string
getnumeric <- function(xstr) as.numeric(unlist(strsplit(xstr, "")))
x <- getnumeric(xstr)
y <- getnumeric(ystr)
#multiply each pair of digits together
mat <- apply(x %o% y, 1, as.character)
#loop over columns, then rows, adding zeroes to end of each number in the matrix to get the correct positioning
ncols <- ncol(mat)
cols <- seq_len(ncols)
for(j in cols)
{
zeroes <- paste(rep("0", ncols-j), collapse="")
mat[,j] <- paste(mat[,j], zeroes, sep="")
}
nrows <- nrow(mat)
rows <- seq_len(nrows)
for(i in rows)
{
zeroes <- paste(rep("0", nrows-i), collapse="")
mat[i,] <- paste(mat[i,], zeroes, sep="")
}
#add zeroes to the start of the each number, so they are all the same length
len <- max(nchar(mat))
strcolumns <- formatC(cbind(as.vector(mat)), width=len)
strcolumns <- gsub(" ", "0", strcolumns)
#line up all the numbers below each other
strmat <- matrix(unlist(strsplit(strcolumns, "")), byrow=TRUE, ncol=len)
#convert to numeric and add them
mat2 <- apply(strmat, 2, as.numeric)
sum1 <- colSums(mat2)
#repeat the process on each of the totals, until each total is a single digit
repeat
{
ntotals <- length(sum1)
totals <- seq_len(ntotals)
for(i in totals)
{
zeroes <- paste(rep("0", ntotals-i), collapse="")
sum1[i] <- paste(sum1[i], zeroes, sep="")
}
len2 <- max(nchar(sum1))
strcolumns2 <- formatC(cbind(as.vector(sum1)), width=len2)
strcolumns2 <- gsub(" ", "0", strcolumns2)
strmat2 <- matrix(unlist(strsplit(strcolumns2, "")), byrow=TRUE, ncol=len2)
mat3 <- apply(strmat2, 2, as.numeric)
sum1 <- colSums(mat3)
if(all(sum1 < 10)) break
}
#Concatenate the digits together
ans <- paste(sum1, collapse="")
ans
}
a <- "18446744073709551616" longmult(a, a)</lang>
"340282366920938463463374607431768211456"
REXX
<lang REXX>/* long multiply */ numeric digits 100 say 2**64*2**64</lang>
Ruby
<lang ruby>def longmult(x,y)
digits = reverse_split_number(x)
result = [0]
j = 0
reverse_split_number(y).each do |m|
c = 0
i = j
digits.each do |d|
v = result[i]
result << 0 if v.zero?
c, v = (v + c + d*m).divmod(10)
result[i] = v
i += 1
end
result[i] += c
j += 1
end
# calculate the answer from the result array of digits
result.reverse.inject(0) {|sum, n| 10*sum + n}
end
def reverse_split_number(m)
digits = [] while m > 0 m, v = m.divmod 10 digits << v end digits
end
n=2**64 printf " %d * %d = %d\n", n, n, n*n printf "longmult(%d, %d) = %d\n", n, n, longmult(n,n)</lang>
18446744073709551616 * 18446744073709551616 = 340282366920938463463374607431768211456 longmult(18446744073709551616, 18446744073709551616) = 340282366920938463463374607431768211456
Scala
This implementation does not rely on an arbitrary precision numeric type. Instead, only single digits are ever multiplied or added, and all partial results are kept as string.
<lang scala>def addNums(x: String, y: String) = {
val padSize = x.length max y.length
val paddedX = "0" * (padSize - x.length) + x
val paddedY = "0" * (padSize - y.length) + y
val (sum, carry) = (paddedX zip paddedY).foldRight(("", 0)) {
case ((dx, dy), (acc, carry)) =>
val sum = dx.asDigit + dy.asDigit + carry
((sum % 10).toString + acc, sum / 10)
}
if (carry != 0) carry.toString + sum else sum
}
def multByDigit(num: String, digit: Int) = {
val (mult, carry) = num.foldRight(("", 0)) {
case (d, (acc, carry)) =>
val mult = d.asDigit * digit + carry
((mult % 10).toString + acc, mult / 10)
}
if (carry != 0) carry.toString + mult else mult
}
def mult(x: String, y: String) =
y.foldLeft("")((acc, digit) => addNums(acc + "0", multByDigit(x, digit.asDigit)))</lang>
Sample:
scala> mult("18446744073709551616", "18446744073709551616")
res25: java.lang.String = 340282366920938463463374607431768211456
Scala 2.8 introduces `scanLeft` and `scanRight` which can be used to simplify this further:
<lang scala>def adjustResult(result: IndexedSeq[Int]) = (
result .map(_ % 10) // remove carry from each digit .tail // drop the seed carry .reverse // put most significant digits on the left .dropWhile(_ == 0) // remove leading zeroes .mkString
)
def addNums(x: String, y: String) = {
val padSize = (x.length max y.length) + 1 // We want to keep a zero to the left, to catch the carry
val paddedX = "0" * (padSize - x.length) + x
val paddedY = "0" * (padSize - y.length) + y
adjustResult((paddedX zip paddedY).scanRight(0) {
case ((dx, dy), last) => dx.asDigit + dy.asDigit + last / 10
})
}
def multByDigit(num: String, digit: Int) = adjustResult(("0"+num).scanRight(0)(_.asDigit * digit + _ / 10))
def mult(x: String, y: String) =
y.foldLeft("")((acc, digit) => addNums(acc + "0", multByDigit(x, digit.asDigit)))
</lang>
Scheme
<lang scheme>(* (expt 2 64) (expt 2 64))</lang>
Seed7
<lang seed7>$ include "seed7_05.s7i";
include "bigint.s7i";
const proc: main is func
begin writeln(2_**64 * 2_**64); end func;</lang>
Output:
340282366920938463463374607431768211456
Slate
<lang slate>(2 raisedTo: 64) * (2 raisedTo: 64).</lang>
Smalltalk
<lang smalltalk>(2 raisedTo: 64) * (2 raisedTo: 64).</lang>
Tcl
Tcl 8.5 supports arbitrary-precision integers, which improves math operations on large integers. It is easy to define our own by following rules for long multiplication; we can then check this against the built-in's result: <lang tcl>package require Tcl 8.5
proc longmult {x y} {
set digits [lreverse [split $x ""]]
set result {0}
set j -2
foreach m [lreverse [split $y ""]] {
set c 0 set i [incr j] foreach d $digits { set v [lindex $result [incr i]] if {$v eq ""} { lappend result 0 set v 0 } regexp (.)(.)$ 0[expr {$v + $c + $d*$m}] -> c v lset result $i $v } lappend result $c
}
# Reconvert digit list into a decimal number
set result [string trimleft [join [lreverse $result] ""] 0]
if {$result == ""} then {return 0} else {return $result}
}
puts [set n [expr {2**64}]] puts [longmult $n $n] puts [expr {$n * $n}]</lang> outputs
18446744073709551616 340282366920938463463374607431768211456 340282366920938463463374607431768211456
Ursala
Natural numbers of unlimited size are a built in type, and arithmetic operations on them are available as library functions. However, since the task calls for explicitly implementing long multiplication, here is an implementation using nothing but language primitives. The numbers are represented as lists of booleans, LSB first. The compiler already knows how to parse and display them in decimal.
<lang Ursala>successor = ~&a^?\1! ~&ah?/~&NfatPRC ~&NNXatPC
sum = ~&B^?a\~&Y@a ~&B?abh/successor@alh2fabt2RC ~&Yabh2Ofabt2RC
product = ~&alrB^& sum@NfalrtPXPRCarh2alPNQX
x = 18446744073709551616
- show+
y = %nP product@iiX x</lang> output:
340282366920938463463374607431768211456