Long multiplication
From Rosetta Code
You are encouraged to solve this task according to the task description, using any language you may know.
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
[edit] Ada
[edit] Using properly range-checked integers
(The source text for these examples can also be found on Bitbucket.)
First we specify the required operations and declare our number type as an array of digits (in base 2^16):
package Long_Multiplication is
type Number (<>) is private;
Zero : constant Number;
One : constant Number;
function Value (Item : in String) return Number;
function Image (Item : in Number) return String;
overriding
function "=" (Left, Right : in Number) return Boolean;
function "+" (Left, Right : in Number) return Number;
function "*" (Left, Right : in Number) return Number;
function Trim (Item : in Number) return Number;
private
Bits : constant := 16;
Base : constant := 2 ** Bits;
type Accumulated_Value is range 0 .. (Base - 1) * Base;
subtype Digit is Accumulated_Value range 0 .. Base - 1;
type Number is array (Natural range <>) of Digit;
for Number'Component_Size use Bits; -- or pragma Pack (Number);
Zero : constant Number := (1 .. 0 => 0);
One : constant Number := (0 => 1);
procedure Divide (Dividend : in Number;
Divisor : in Digit;
Result : out Number;
Remainder : out Digit);
end Long_Multiplication;
Some of the operations declared above are useful helper operations for the conversion of numbers to and from base 10 digit strings.
Then we implement the operations:
package body Long_Multiplication is
function Value (Item : in String) return Number is
subtype Base_Ten_Digit is Digit range 0 .. 9;
Ten : constant Number := (0 => 10);
begin
case Item'Length is
when 0 =>
raise Constraint_Error;
when 1 =>
return (0 => Base_Ten_Digit'Value (Item));
when others =>
return (0 => Base_Ten_Digit'Value (Item (Item'Last .. Item'Last)))
+ Ten * Value (Item (Item'First .. Item'Last - 1));
end case;
end Value;
function Image (Item : in Number) return String is
Base_Ten : constant array (Digit range 0 .. 9) of String (1 .. 1) :=
("0", "1", "2", "3", "4", "5", "6", "7", "8", "9");
Result : Number (0 .. Item'Last);
Remainder : Digit;
begin
if Item = Zero then
return "0";
else
Divide (Dividend => Item,
Divisor => 10,
Result => Result,
Remainder => Remainder);
if Result = Zero then
return Base_Ten (Remainder);
else
return Image (Trim (Result)) & Base_Ten (Remainder);
end if;
end if;
end Image;
overriding
function "=" (Left, Right : in Number) return Boolean is
begin
for Position in Integer'Min (Left'First, Right'First) ..
Integer'Max (Left'Last, Right'Last) loop
if Position in Left'Range and Position in Right'Range then
if Left (Position) /= Right (Position) then
return False;
end if;
elsif Position in Left'Range then
if Left (Position) /= 0 then
return False;
end if;
elsif Position in Right'Range then
if Right (Position) /= 0 then
return False;
end if;
else
raise Program_Error;
end if;
end loop;
return True;
end "=";
function "+" (Left, Right : in Number) return Number is
Result : Number (Integer'Min (Left'First, Right'First) ..
Integer'Max (Left'Last , Right'Last) + 1);
Accumulator : Accumulated_Value := 0;
Used : Integer := Integer'First;
begin
for Position in Result'Range loop
if Position in Left'Range then
Accumulator := Accumulator + Left (Position);
end if;
if Position in Right'Range then
Accumulator := Accumulator + Right (Position);
end if;
Result (Position) := Accumulator mod Base;
Accumulator := Accumulator / Base;
if Result (Position) /= 0 then
Used := Position;
end if;
end loop;
if Accumulator = 0 then
return Result (Result'First .. Used);
else
raise Constraint_Error;
end if;
end "+";
function "*" (Left, Right : in Number) return Number is
Accumulator : Accumulated_Value;
Result : Number (Left'First + Right'First ..
Left'Last + Right'Last + 1) := (others => 0);
Used : Integer := Integer'First;
begin
for L in Left'Range loop
for R in Right'Range loop
Accumulator := Left (L) * Right (R);
for Position in L + R .. Result'Last loop
exit when Accumulator = 0;
Accumulator := Accumulator + Result (Position);
Result (Position) := Accumulator mod Base;
Accumulator := Accumulator / Base;
Used := Position;
end loop;
end loop;
end loop;
return Result (Result'First .. Used);
end "*";
procedure Divide (Dividend : in Number;
Divisor : in Digit;
Result : out Number;
Remainder : out Digit) is
Accumulator : Accumulated_Value := 0;
begin
Result := (others => 0);
for Position in reverse Dividend'Range loop
Accumulator := Accumulator * Base + Dividend (Position);
Result (Position) := Accumulator / Divisor;
Accumulator := Accumulator mod Divisor;
end loop;
Remainder := Accumulator;
end Divide;
function Trim (Item : in Number) return Number is
begin
for Position in reverse Item'Range loop
if Item (Position) /= 0 then
return Item (Item'First .. Position);
end if;
end loop;
return Zero;
end Trim;
end Long_Multiplication;
And finally we have the requested test application:
with Ada.Text_IO;
with Long_Multiplication;
procedure Test_Long_Multiplication is
use Ada.Text_IO, Long_Multiplication;
N : Number := Value ("18446744073709551616");
M : Number := N * N;
begin
Put_Line (Image (N) & " * " & Image (N) & " = " & Image (M));
end Test_Long_Multiplication;
- Output:
18446744073709551616 * 18446744073709551616 = 340282366920938463463374607431768211456
[edit] Using modular types
The following implementation uses representation of a long number by an array of 32-bit elements:
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 "*";
The task requires conversion into decimal base. For this we also need division to short number with a remainder. Here it is:
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;
With the above the test program:
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;
Sample output:
340282366920938463463374607431768211456
[edit] 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.
[edit] 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.
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))
)
Output:
The approximate golden ratios, width: +618033988749894848204586834365638117720309179805762862135448622705261
length: +1618033988749894848204586834365638117720309179805762862135448622705261
product is exactly: +1000000000000000000000000000000000000000000000000000000000000000000001201173450350400438606015942314498798603569682901026716145698077078121
2 ** 64 * -(2 ** 64) = -340282366920938463463374607431768211456
[edit] Implementation example
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
);
################################################################
# 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
);
# 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))
)
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
[edit] Other libraries or implementation specific extensions
As of February 2009 no open source libraries to do this task have been located.
[edit] AutoHotkey
ahk discussion
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)
}
[edit] 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 "===="
}
}
outputs:
2^64 * 2^64 = 340282366920938463463374607431768211456 2^64 * 2^64 = 340282366920938463463374607431768211456
[edit] BASIC
[edit] Version 1
'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
[edit] Version 2
'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"
[edit] BBC BASIC
Library method:
INSTALL @lib$+"BB4WMAPMLIB"
MAPM_DllPath$ = @lib$+"BB4WMAPM.DLL"
PROCMAPM_Init
twoto64$ = "18446744073709551616"
PRINT "2^64 * 2^64 = " ; FNMAPM_Multiply(twoto64$, twoto64$)
Explicit method:
twoto64$ = "18446744073709551616"
PRINT "2^64 * 2^64 = " ; FNlongmult(twoto64$, twoto64$)
END
DEF FNlongmult(num1$, num2$)
LOCAL C%, I%, J%, S%, num1&(), num2&(), num3&()
S% = LEN(num1$)+LEN(num2$)
DIM num1&(S%), num2&(S%), num3&(S%)
IF LEN(num1$) > LEN(num2$) SWAP num1$,num2$
$$^num1&(1) = num1$
num1&() AND= 15
FOR I% = LEN(num1$) TO 1 STEP -1
$$^num2&(I%) = num2$
num2&() AND= 15
num3&() += num2&() * num1&(I%)
IF I% MOD 3 = 1 THEN
C% = 0
FOR J% = S%-1 TO I%-1 STEP -1
C% += num3&(J%)
num3&(J%) = C% MOD 10
C% DIV= 10
NEXT
ENDIF
NEXT I%
num3&() += &30
num3&(S%) = 0
IF num3&(0) = &30 THEN = $$^num3&(1)
= $$^num3&(0)
[edit] C
Doing it as if by hand.
#include <stdio.h>output
#include <string.h>
/* c = a * b. Caller is responsible for memory.
c must not be the same as either a or b. */
void longmulti(const char *a, const char *b, char *c)
{
int i = 0, j = 0, k = 0, n, carry;
int la, lb;
/* either is zero, return "0" */
if (!strcmp(a, "0") || !strcmp(b, "0")) {
c[0] = '0', c[1] = '\0';
return;
}
/* see if either a or b is negative */
if (a[0] == '-') { i = 1; k = !k; }
if (b[0] == '-') { j = 1; k = !k; }
/* if yes, prepend minus sign if needed and skip the sign */
if (i || j) {
if (k) c[0] = '-';
longmulti(a + i, b + j, c + k);
return;
}
la = strlen(a);
lb = strlen(b);
memset(c, '0', la + lb);
c[la + lb] = '\0';
# define I(a) (a - '0')
for (i = la - 1; i >= 0; i--) {
for (j = lb - 1, k = i + j + 1, carry = 0; j >= 0; j--, k--) {
n = I(a[i]) * I(b[j]) + I(c[k]) + carry;
carry = n / 10;
c[k] = (n % 10) + '0';
}
c[k] += carry;
}
# undef I
if (c[0] == '0') memmove(c, c + 1, la + lb);
return;
}
int main()
{
char c[1024];
longmulti("-18446744073709551616", "-18446744073709551616", c);
printf("%s\n", c);
return 0;
}
340282366920938463463374607431768211456
[edit] C++
[edit] Version 1
#include <iostream>
#include <sstream>
//--------------------------------------------------------------------------------------------------
typedef long long bigInt;
//--------------------------------------------------------------------------------------------------
using namespace std;
//--------------------------------------------------------------------------------------------------
class number
{
public:
number() { s = "0"; neg = false; }
number( bigInt a ) { set( a ); }
number( string a ) { set( a ); }
void set( bigInt a ) { neg = false; if( a < 0 ) { a = -a; neg = true; } ostringstream o; o << a; s = o.str(); clearStr(); }
void set( string a ) { neg = false; s = a; if( s.length() > 1 && s[0] == '-' ) { neg = true; } clearStr(); }
number operator * ( const number& b ) { return this->mul( b ); }
number& operator *= ( const number& b ) { *this = *this * b; return *this; }
number& operator = ( const number& b ) { s = b.s; return *this; }
friend ostream& operator << ( ostream& out, const number& a ) { if( a.neg ) out << "-"; out << a.s; return out; }
friend istream& operator >> ( istream& in, number& a ){ string b; in >> b; a.set( b ); return in; }
private:
number mul( const number& b )
{
number a; bool neg = false;
string r, bs = b.s; r.resize( 2 * max( b.s.length(), s.length() ), '0' );
int xx, ss, rr, t, c, stp = 0;
string::reverse_iterator xi = bs.rbegin(), si, ri;
for( ; xi != bs.rend(); xi++ )
{
c = 0; ri = r.rbegin() + stp;
for( si = s.rbegin(); si != s.rend(); si++ )
{
xx = ( *xi ) - 48; ss = ( *si ) - 48; rr = ( *ri ) - 48;
ss = ss * xx + rr + c; t = ss % 10; c = ( ss - t ) / 10;
( *ri++ ) = t + 48;
}
if( c > 0 ) ( *ri ) = c + 48;
stp++;
}
trimLeft( r ); t = b.neg ? 1 : 0; t += neg ? 1 : 0;
if( t & 1 ) a.s = "-" + r;
else a.s = r;
return a;
}
void trimLeft( string& r )
{
if( r.length() < 2 ) return;
for( string::iterator x = r.begin(); x != ( r.end() - 1 ); )
{
if( ( *x ) != '0' ) return;
x = r.erase( x );
}
}
void clearStr()
{
for( string::iterator x = s.begin(); x != s.end(); )
{
if( ( *x ) < '0' || ( *x ) > '9' ) x = s.erase( x );
else x++;
}
}
string s;
bool neg;
};
//--------------------------------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
number a, b;
a.set( "18446744073709551616" ); b.set( "18446744073709551616" );
cout << a * b << endl << endl;
cout << "Factor 1 = "; cin >> a;
cout << "Factor 2 = "; cin >> b;
cout << "Product: = " << a * b << endl << endl;
return system( "pause" );
}
//--------------------------------------------------------------------------------------------------
- Output:
340282366920938463463374607431768211456 Factor 1 = 9876548974569852365985574874787454878778975948 Factor 2 = 8954564845421878741168741154541897945138974567 Product: = 88440198241770705041777453160463400993104404280916080859287340887463980926235972531076714516
[edit] Version 2
using the Class Library for Numbers , CLN, and g++
#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 ;
}
Output: The result of 2^64 * 2^64 is 340282366920938463463374607431768211456 !
[edit] C#
System.Numerics.BigInteger was added with C# 4.
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);
}
}
Output:
340282366920938463463374607431768211456
[edit] CoffeeScript
# This very limited BCD-based collection of functions
# allows for long multiplication. It works for positive
# numbers only. The assumed data structure is as follows:
# BcdInteger.from_integer(4321) == [1, 2, 3, 4]
BcdInteger =
from_string: (s) ->
arr = []
for c in s
arr.unshift parseInt(c)
arr
from_integer: (n) ->
result = []
while n > 0
result.push n % 10
n = Math.floor n / 10
result
to_string: (arr) ->
s = ''
for elem in arr
s = elem.toString() + s
s
sum: (arr1, arr2) ->
if arr1.length < arr2.length
return BcdInteger.sum(arr2, arr1)
carry = 0
result= []
for d1, pos in arr1
d = d1 + (arr2[pos] || 0) + carry
result.push d % 10
carry = Math.floor d / 10
if carry
result.push 1
result
multiply_by_power_of_ten: (arr, power_of_ten) ->
result = (0 for i in [0...power_of_ten])
result.concat arr
product_by_integer: (arr, n) ->
result = []
for digit, i in arr
prod = BcdInteger.from_integer n * digit
prod = BcdInteger.multiply_by_power_of_ten prod, i
result = BcdInteger.sum result, prod
result
product: (arr1, arr2) ->
result = []
for digit, i in arr1
prod = BcdInteger.product_by_integer arr2, digit
prod = BcdInteger.multiply_by_power_of_ten prod, i
result = BcdInteger.sum result, prod
result
x = BcdInteger.from_integer 1
for i in [1..64]
x = BcdInteger.product_by_integer x, 2
console.log BcdInteger.to_string x # 18446744073709551616
square = BcdInteger.product x, x
console.log BcdInteger.to_string square # 340282366920938463463374607431768211456
[edit] Common 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)))))
> (long-multiply (expt 2 64) (expt 2 64)) 340282366920938463463374607431768211456
[edit] D
Using the standard library:
void main() {
import std.stdio, std.bigint;
writeln(2.BigInt ^^ 64 * 2.BigInt ^^ 64);
}
- Output:
340282366920938463463374607431768211456
Long multiplication, same output:
import std.stdio, std.algorithm, std.range;
auto longMult(in string x1, in string x2) /*pure nothrow*/ {
auto digits1 = x1.retro.map!q{a - '0'};
const digits2 = x2.retro.map!q{a - '0'}.array;
uint[] res;
foreach (immutable i, immutable d1; uint.max.iota.zip(digits1)) {
foreach (immutable j, immutable d2; digits2) {
immutable k = i + j;
if (res.length <= k)
res.length++;
res[k] += d1 * d2;
if (res[k] > 9) {
if (res.length <= k + 1)
res.length++;
res[k + 1] = res[k] / 10 + res[k + 1];
res[k] -= res[k] / 10 * 10;
}
}
}
return res.retro.map!q{ cast(char)(a + '0') };
}
void main() {
immutable two64 = "18446744073709551616";
longMult(two64, two64).writeln;
}
[edit] Dc
Since Dc has arbitrary precision built-in, the task is no different than a normal multiplication:
2 64^ 2 64^ *p
[edit] 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)})
Output:
a is 4294967296 a*a is 18446744073709551616 a*a*a*a is 340282366920938463488374607424768211456
[edit] F#
> let X = 2I ** 64 * 2I ** 64 ;;
val X : System.Numerics.BigInteger = 340282366920938463463374607431768211456
[edit] Factor
USING: kernel math sequences ;
: longmult-seq ( xs ys -- zs )
[ * ] cartesian-map
dup length iota [ 0 <repetition> ] map
[ prepend ] 2map
[ ] [ [ 0 suffix ] dip [ + ] 2map ] map-reduce ;
: integer->digits ( x -- xs ) { } swap [ dup 0 > ] [ 10 /mod swap [ prefix ] dip ] while drop ;
: digits->integer ( xs -- x ) 0 [ swap 10 * + ] reduce ;
: longmult ( x y -- z ) [ integer->digits ] bi@ longmult-seq digits->integer ;
( scratchpad ) 2 64 ^ dup longmult .
340282366920938463463374607431768211456
( scratchpad ) 2 64 ^ dup * .
340282366920938463463374607431768211456
[edit] 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
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
[edit] Go
// Long multiplication per WP article referenced by task description.
// That is, multiplicand is multiplied by single digits of multiplier
// to form intermediate results. Intermediate results are accumulated
// for the product. Used here is the abacus method mentioned by the
// article, of summing intermediate results as they are produced,
// rather than all at once at the end.
//
// Limitations: Negative numbers not supported, superfluous leading zeros
// not generally removed.
package main
import "fmt"
// argument validation
func d(b byte) byte {
if b < '0' || b > '9' {
panic("digit 0-9 expected")
}
return b - '0'
}
// add two numbers as strings
func add(x, y string) string {
if len(y) > len(x) {
x, y = y, x
}
b := make([]byte, len(x)+1)
var c byte
for i := 1; i <= len(x); i++ {
if i <= len(y) {
c += d(y[len(y)-i])
}
s := d(x[len(x)-i]) + c
c = s / 10
b[len(b)-i] = (s % 10) + '0'
}
if c == 0 {
return string(b[1:])
}
b[0] = c + '0'
return string(b)
}
// multipy a number by a single digit
func mulDigit(x string, y byte) string {
if y == '0' {
return "0"
}
y = d(y)
b := make([]byte, len(x)+1)
var c byte
for i := 1; i <= len(x); i++ {
s := d(x[len(x)-i])*y + c
c = s / 10
b[len(b)-i] = (s % 10) + '0'
}
if c == 0 {
return string(b[1:])
}
b[0] = c + '0'
return string(b)
}
// multiply two numbers as strings
func mul(x, y string) string {
result := mulDigit(x, y[len(y)-1])
for i, zeros := 2, ""; i <= len(y); i++ {
zeros += "0"
result = add(result, mulDigit(x, y[len(y)-i])+zeros)
}
return result
}
// requested output
const n = "18446744073709551616"
func main() {
fmt.Println(mul(n, n))
}
Output:
340282366920938463463374607431768211456
[edit] Groovy
println 2**64 * 2**64
340282366920938463463374607431768211456
[edit] Haskell
import Data.List (transpose)
import Data.Char (digitToInt)
import Data.List (inits)
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
main = print $ (2 ^ 64) `longmult` (2 ^ 64)
- Output:
340282366920938463463374607431768211456
[edit] Icon and Unicon
Icon and Unicon support large integers.
procedure main()
write(2^64*2^64)
end
[edit] J
Solution:
digits =: ,.&.":
polymult =: +//.@(*/)
buildDecimal=: 10x&#.
longmult=: buildDecimal@polymult&digits
Example:
longmult~ 2x^64
340282366920938463463374607431768211456
Alternatives:
longmult could have been defined concisely:
longmult=: 10x&#.@(+//.@(*/)&(,.&.":))
Or, of course, the task may be accomplished without the verb definitions:
10x&#.@(+//.@(*/)&(,.&.":))~2x^64
340282366920938463463374607431768211456
Or using the code (+ 10x&*)/@|. instead of #.:
(+ 10x&*)/@|.@(+//.@(*/)&(,.&.":))~2x^64
340282366920938463463374607431768211456
Or you could use the built-in language support for arbitrary precision multiplication:
(2x^64)*(2x^64)
340282366920938463463374607431768211456
Explaining the component verbs:
-
digitstranslates a number to a corresponding list of digits;
,.&.": 123
1 2 3
-
polymult(multiplies polynomials): ref. [1]
1 2 3 (+//.@(*/)) 1 2 3
1 4 10 12 9
-
buildDecimal(translates a list of decimal digits to the corresponding extended precision number):
(+ 10x&*)/|. 1 4 10 12 9
15129
[edit] Java
[edit] Decimal version
This version of the code keeps the data in base ten. By doing this, we can avoid converting the whole number to binary and we can keep things simple, but the runtime will be suboptimal.
public class LongMult {
private static byte[] stringToDigits(String num) {
byte[] result = new byte[num.length()];
for (int i = 0; i < num.length(); i++) {
char c = num.charAt(i);
if (c < '0' || c > '9') {
throw new IllegalArgumentException("Invalid digit " + c
+ " found at position " + i);
}
result[num.length() - 1 - i] = (byte) (c - '0');
}
return result;
}
public static String longMult(String num1, String num2) {
byte[] left = stringToDigits(num1);
byte[] right = stringToDigits(num2);
byte[] result = new byte[left.length + right.length];
for (int rightPos = 0; rightPos < right.length; rightPos++) {
byte rightDigit = right[rightPos];
byte temp = 0;
for (int leftPos = 0; leftPos < left.length; leftPos++) {
temp += result[leftPos + rightPos];
temp += rightDigit * left[leftPos];
result[leftPos + rightPos] = (byte) (temp % 10);
temp /= 10;
}
int destPos = rightPos + left.length;
while (temp != 0) {
temp += result[destPos] & 0xFFFFFFFFL;
result[destPos] = (byte) (temp % 10);
temp /= 10;
destPos++;
}
}
StringBuilder stringResultBuilder = new StringBuilder(result.length);
for (int i = result.length - 1; i >= 0; i--) {
byte digit = result[i];
if (digit != 0 || stringResultBuilder.length() > 0) {
stringResultBuilder.append((char) (digit + '0'));
}
}
return stringResultBuilder.toString();
}
public static void main(String[] args) {
System.out.println(longMult("18446744073709551616",
"18446744073709551616"));
}
}
[edit] Binary version
This version tries to be as efficient as possible, so it converts numbers into binary before doing any calculations. The complexity is higher because of the need to convert to and from base ten, which requires the implementation of some additional arithmetic operations beyond long multiplication itself.
import java.util.Arrays;
public class LongMultBinary {
/**
* A very basic arbitrary-precision integer class. It only handles
* non-negative numbers and doesn't implement any arithmetic not necessary
* for the task at hand.
*/
public static class MyLongNum implements Cloneable {
/*
* The actual bits of the integer, with the least significant place
* first. The biggest native integer type of Java is the 64-bit long,
* but since we need to be able to store the result of two digits
* multiplied, we have to use the second biggest native type, the 32-bit
* int. All numeric types are signed in Java, but we don't want to waste
* the sign bit, so we need to take extra care while doing arithmetic to
* ensure unsigned semantics.
*/
private int[] digits;
/*
* The number of digits actually used in the digits array. Since arrays
* cannot be resized in Java, we are better off remembering the logical
* size ourselves, instead of reallocating and copying every time we need to shrink.
*/
private int digitsUsed;
@Override
public MyLongNum clone() {
try {
MyLongNum clone = (MyLongNum) super.clone();
clone.digits = clone.digits.clone();
return clone;
} catch (CloneNotSupportedException e) {
throw new Error("Object.clone() threw exception", e);
}
}
private void resize(int newLength) {
if (digits.length < newLength) {
digits = Arrays.copyOf(digits, newLength);
}
}
private void adjustDigitsUsed() {
while (digitsUsed > 0 && digits[digitsUsed - 1] == 0) {
digitsUsed--;
}
}
/**
* "Short" multiplication by one digit. Used to convert strings to long numbers.
*/
public void multiply(int multiplier) {
if (multiplier < 0) {
throw new IllegalArgumentException(
"Signed arithmetic isn't supported");
}
resize(digitsUsed + 1);
long temp = 0;
for (int i = 0; i < digitsUsed; i++) {
temp += (digits[i] & 0xFFFFFFFFL) * multiplier;
digits[i] = (int) temp; // store the low 32 bits
temp >>>= 32;
}
digits[digitsUsed] = (int) temp;
digitsUsed++;
adjustDigitsUsed();
}
/**
* "Short" addition (adding a one-digit number). Used to convert strings to long numbers.
*/
public void add(int addend) {
if (addend < 0) {
throw new IllegalArgumentException(
"Signed arithmetic isn't supported");
}
long temp = addend;
for (int i = 0; i < digitsUsed && temp != 0; i++) {
temp += (digits[i] & 0xFFFFFFFFL);
digits[i] = (int) temp; // store the low 32 bits
temp >>>= 32;
}
if (temp != 0) {
resize(digitsUsed + 1);
digits[digitsUsed] = (int) temp;
digitsUsed++;
}
}
/**
* "Short" division (dividing by a one-digit number). Used to convert numbers to strings.
* @param divisor The digit to divide by.
* @return The remainder of the division.
*/
public int divide(int divisor) {
if (divisor < 0) {
throw new IllegalArgumentException(
"Signed arithmetic isn't supported");
}
int remainder = 0;
for (int i = digitsUsed - 1; i >= 0; i--) {
long twoDigits = (((long) remainder << 32) | (digits[i] & 0xFFFFFFFFL));
remainder = (int) (twoDigits % divisor);
digits[i] = (int) (twoDigits / divisor);
}
adjustDigitsUsed();
return remainder;
}
public MyLongNum(String value) {
// each of our 32-bit digits can store at least 9 decimal digit's worth
this.digits = new int[value.length() / 9 + 1];
this.digitsUsed = 0;
// To lower the number of bignum operations, handle nine digits at a time.
for (int i = 0; i < value.length(); i+=9) {
String chunk = value.substring(i, Math.min(i+9, value.length()));
int multiplier = 1;
int addend = 0;
for (int j=0; j<chunk.length(); j++) {
char c = chunk.charAt(j);
if (c < '0' || c > '9') {
throw new IllegalArgumentException("Invalid digit " + c
+ " found in input");
}
multiplier *= 10;
addend *= 10;
addend += c - '0';
}
multiply(multiplier);
add(addend);
}
}
@Override
public String toString() {
if (digitsUsed == 0) {
return "0";
}
MyLongNum dummy = this.clone();
StringBuilder resultBuilder = new StringBuilder(digitsUsed * 9);
while (dummy.digitsUsed > 0) {
// To limit the number of bignum divisions, handle nine digits at a time.
int decimalDigits = dummy.divide(1000000000);
for (int i=0; i<9; i++) {
resultBuilder.append((char) (decimalDigits % 10 + '0'));
decimalDigits /= 10;
}
}
// Trim any leading zeros we may have created.
while (resultBuilder.charAt(resultBuilder.length()-1) == '0') {
resultBuilder.deleteCharAt(resultBuilder.length()-1);
}
return resultBuilder.reverse().toString();
}
/**
* Long multiplication.
*/
public void multiply(MyLongNum multiplier) {
MyLongNum left, right;
// Make sure the shorter number is on the right-hand side to make things a bit more efficient.
if (this.digitsUsed > multiplier.digitsUsed) {
left = this;
right = multiplier;
} else {
left = multiplier;
right = this;
}
int[] newDigits = new int[left.digitsUsed + right.digitsUsed];
for (int rightPos = 0; rightPos < right.digitsUsed; rightPos++) {
long rightDigit = right.digits[rightPos] & 0xFFFFFFFFL;
long temp = 0;
for (int leftPos = 0; leftPos < left.digitsUsed; leftPos++) {
temp += (newDigits[leftPos + rightPos] & 0xFFFFFFFFL);
temp += rightDigit * (left.digits[leftPos] & 0xFFFFFFFFL);
newDigits[leftPos + rightPos] = (int) temp;
temp >>>= 32;
}
// Roll forward any carry we may have.
int destPos = rightPos + digitsUsed;
while (temp != 0) {
temp += (newDigits[destPos] & 0xFFFFFFFFL);
newDigits[destPos] = (int) temp;
temp >>>= 32;
destPos++;
}
}
this.digits = newDigits;
this.digitsUsed = newDigits.length;
adjustDigitsUsed();
}
}
public static void main(String[] args) {
MyLongNum one = new MyLongNum("18446744073709551616");
MyLongNum two = one.clone();
one.multiply(two);
System.out.println(one);
}
}
[edit] 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("");
}
[edit] Liberty BASIC
'[RC] long multiplication
'now, count 2^64
print "2^64"
a$="1"
for i = 1 to 64
a$ = multByD$(a$, 2)
next
print a$
print "(check with native LB)"
print 2^64
print "(looks OK)"
'now let's do b$*a$ stuff
print "2^64*2^64"
print longMult$(a$, a$)
print "(check with native LB)"
print 2^64*2^64
print "(looks OK)"
end
'---------------------------------------
function longMult$(a$, b$)
signA = 1
if left$(a$,1) = "-" then a$ = mid$(a$,2): signA = -1
signB = 1
if left$(b$,1) = "-" then b$ = mid$(b$,2): signB = -1
c$ = ""
t$ = ""
shift$ = ""
for i = len(a$) to 1 step -1
d = val(mid$(a$,i,1))
t$ = multByD$(b$, d)
c$ = addLong$(c$, t$+shift$)
shift$ = shift$ +"0"
'print d, t$, c$
next
if signA*signB<0 then c$ = "-" + c$
'print c$
longMult$ = c$
end function
function multByD$(a$, d)
'multiply a$ by digit d
c$ = ""
carry = 0
for i = len(a$) to 1 step -1
a = val(mid$(a$,i,1))
c = a*d+carry
carry = int(c/10)
c = c mod 10
'print a, c
c$ = str$(c)+c$
next
if carry>0 then c$ = str$(carry)+c$
'print c$
multByD$ = c$
end function
function addLong$(a$, b$)
'add a$ + b$, for now only positive
l = max(len(a$), len(b$))
a$=pad$(a$,l)
b$=pad$(b$,l)
c$ = "" 'result
carry = 0
for i = l to 1 step -1
a = val(mid$(a$,i,1))
b = val(mid$(b$,i,1))
c = a+b+carry
carry = int(c/10)
c = c mod 10
'print a, b, c
c$ = str$(c)+c$
next
if carry>0 then c$ = str$(carry)+c$
'print c$
addLong$ = c$
end function
function pad$(a$,n) 'pad from right with 0 to length n
pad$ = a$
while len(pad$)<n
pad$ = "0"+pad$
wend
end function
[edit] Mathematica
We define the long multiplication function:
LongMultiplication[a_,b_]:=Module[{d1,d2},
d1=IntegerDigits[a]//Reverse;
d2=IntegerDigits[b]//Reverse;
Sum[d1[[i]]d2[[j]]*10^(i+j-2),{i,1,Length[d1]},{j,1,Length[d2]}]
]
Example:
n1 = 2^64;
n2 = 2^64;
LongMultiplication[n1, n2]
gives back:
340282366920938463463374607431768211456
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:
n1=2^8000;
n2=2^8000;
Timing[LongMultiplication[n1,n2]][[1]]
Timing[n1 n2][[1]]
Floor[%%/%]
gives back:
72.9686
7.*10^-6
10424088
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.
[edit] PARI/GP
long(a,b)={
a=eval(Vec(a));
b=eval(Vec(b));
my(c=vector(#a+#b),carry=0);
for(i=1,#a,
for(j=1,#b,
c[i+j]+=a[i]*b[j]
)
);
forstep(i=#c,1,-1,
c[i] += carry;
carry = c[i] \ 10;
c[i] = c[i] % 10
);
for(i=1,#c,
if(c[i], return(concat(apply(s->Str(s),vector(#c+1-i,j,c[i+j-1])))))
);
"0"
};
long("18446744073709551616","18446744073709551616")
Output:
%1 = "340282366920938463463374607431768211456"
[edit] 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";
[edit] 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.
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;
In any case, integers are specced to be arbitrarily large, which at least two implementations (Pugs and Niecza) support already:
niecza> 18446744073709551616 * 18446744073709551616 340282366920938463463374607431768211456
[edit] PHP
<?php
$factor = bcpow(2, 64);
$product = bcmul($factor, $factor);
echo "2^64 * 2^64 is " . $product;
?>
Output:
2^64 * 2^64 is 340282366920938463463374607431768211456
[edit] PicoLisp
: (* (** 2 64) (** 2 64))
-> 340282366920938463463374607431768211456
[edit] 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;
Calling sequence:
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);
Final output:
18446744073709551616 340282366920938463463374607431768211456
[edit] 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"
[edit] Prolog
Arbitrary precision arithmetic is native in most Prolog implementations.
?- X is 2**64 * 2**64.
X = 340282366920938463463374607431768211456.
[edit] PureBasic
[edit] Explicit Implementation
Structure decDigitFmt ;decimal digit format
Array Digit.b(0) ;contains each digit of number, right-most digit is index 0
digitCount.i ;zero based
sign.i ; {x < 0} = -1, {x = 0} = 0, {x > 0} = 1
EndStructure
Global zero_decDigitFmt.decDigitFmt ;represents zero in the decimal digit format
;converts string representation of integer into the digit format, number can include signus but no imbedded spaces
Procedure stringToDecDigitFmt(numString.s, *x.decDigitFmt)
Protected *c.Character, digitIdx, digitCount
If numString And *x
*c.Character = @numString
Repeat
Select *c\c
Case '0' To '9', '-', '+'
*c + SizeOf(Character)
Default
numString = Left(numString, *c - @numString)
Break
EndSelect
ForEver
*c = @numString
Select *c\c
Case '-'
*x\sign = -1
*c + SizeOf(Character)
Case '+'
*x\sign = 1
*c + SizeOf(Character)
Case '0' To '9'
*x\sign = 1
EndSelect
numString = LTrim(PeekS(*c), "0") ;remove leading zeroes
If numString = "" ;is true if equal to zero or if only a signus is present
CopyStructure(@zero_decDigitFmt, *x, decDigitFmt)
ProcedureReturn
EndIf
*c = @numString
digitCount = Len(PeekS(*c)) - 1
Dim *x\Digit(digitCount)
*x\digitCount = digitCount
digitIdx = 0
While *c\c
If *c\c >= '0' And *c\c <= '9'
*x\Digit(digitCount - digitIdx) = *c\c - '0'
digitIdx + 1
*c + SizeOf(Character)
Else
Break
EndIf
Wend
EndIf
EndProcedure
;converts digit format representation of integer into string representation
Procedure.s decDigitFmtToString(*x.decDigitFmt)
Protected i, number.s
If *x
If *x\sign = 0
number = "0"
Else
For i = *x\digitCount To 0 Step -1
number + Str(*x\Digit(i))
Next
number = LTrim(number, "0")
If *x\sign = -1
number = "-" + number
EndIf
EndIf
EndIf
ProcedureReturn number
EndProcedure
;handles only positive numbers and zero, negative numbers left as an exercise for the reader ;)
Procedure add_decDigitFmt(*a.decDigitFmt, *b.decDigitFmt, *sum.decDigitFmt, digitPos = 0) ;*sum contains the result of (*a ) * 10^digitPos + (*b)
Protected carry, i, newDigitCount, workingSum, a_dup.decDigitFmt
If *a And *b And *sum
If *a = *sum: CopyStructure(*a, @a_dup, decDigitFmt): *a = @a_dup: EndIf ;handle special case of *sum + *b = *sum
If *b <> *sum: CopyStructure(*b, *sum, decDigitFmt): EndIf ;handle general case of *a + *b = *sum and special case of *a + *sum = *sum
;calculate number of digits needed for sum and resize array of digits if necessary
newDigitCount = *a\digitCount + digitPos
If newDigitCount >= *sum\digitCount
If *sum\digitCount = newDigitCount And *sum\Digit(*sum\digitCount) <> 0
newDigitCount + 1
EndIf
If *sum\digitCount <> newDigitCount
*sum\digitCount = newDigitCount
Redim *sum\Digit(*sum\digitCount)
EndIf
EndIf
i = 0
Repeat
If i <= *a\digitCount
workingSum = *a\Digit(i) + *sum\Digit(digitPos) + carry
Else
workingSum = *sum\Digit(digitPos) + carry
EndIf
If workingSum > 9
carry = 1
workingSum - 10
Else
carry = 0
EndIf
*sum\Digit(digitPos) = workingSum
digitPos + 1
i + 1
Until i > *a\digitCount And carry = 0
If *a\sign <> 0 Or *sum\sign <> 0
*sum\sign = 1 ;only handle positive numbers and zero for now
EndIf
EndIf
EndProcedure
Procedure multiply_decDigitFmt(*a.decDigitFmt, *b.decDigitFmt, *product.decDigitFmt) ;*product contains the result of (*a) * (*b)
Protected i, digitPos, productSignus
Protected Dim multTable.decDigitFmt(9)
Protected NewList digitProduct.decDigitFmt()
If *a And *b And *product
If *a\sign = 0 Or *b\sign = 0
CopyStructure(zero_decDigitFmt, *product, decDigitFmt)
ProcedureReturn
EndIf
If *b\digitCount > *a\digitCount: Swap *a, *b: EndIf
;build multiplication table
CopyStructure(*a, @multTable(1), decDigitFmt): multTable(1)\sign = 1 ;always positive
For i = 2 To 9
add_decDigitFmt(*a, multTable(i - 1), multTable(i))
Next
;collect individual digit products for later summation; these could also be added as we go along
For i = 0 To *b\digitCount
AddElement(digitProduct())
digitProduct() = multTable(*b\Digit(i))
Next
;determine sign of product
If *a\sign <> *b\sign
productSignus = -1
Else
productSignus = 1
EndIf
digitPos = 0
CopyStructure(zero_decDigitFmt, *product, decDigitFmt)
ForEach digitProduct()
add_decDigitFmt(digitProduct(), *product, *product, digitPos)
digitPos + 1
Next
*product\sign = productSignus ;set sign of product
EndIf
EndProcedure
;handles only positive integer exponents or an exponent of zero, does not raise an error for 0^0
Procedure exponent_decDigitFmt(*a.decDigitFmt, exponent, *product.decDigitFmt)
Protected i, a_dup.decDigitFmt
If *a And *product And exponent >= 0
If *a = *product: CopyStructure(*a, @a_dup, decDigitFmt): *a = @a_dup: EndIf
stringToDecDigitFmt("1", *product)
For i = 1 To exponent: multiply_decDigitFmt(*product, *a, *product): Next
EndIf
EndProcedure
If OpenConsole()
Define a.decDigitFmt, product.decDigitFmt
stringToDecDigitFmt("2", a)
exponent_decDigitFmt(a, 64, a) ;2^64
multiply_decDigitFmt(a, a, product)
PrintN("The result of 2^64 * 2^64 is " + decDigitFmtToString(product))
Print(#crlf$ + #crlf$ + "Press ENTER to exit"): Input()
CloseConsole()
EndIf
Output:
The result of 2^64 * 2^64 is 340282366920938463463374607431768211456
[edit] Library Method
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.
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))
Outputs
2^64*2^64 = 340282366920938463463374607431768211456
[edit] Python
(Note that Python comes with arbitrary length integers).
#!/usr/bin/env python
print 2**64*2**64
#!/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)
Shorter version:
#!/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)
[edit] R
[edit] Using GMP
library(gmp)
a <- as.bigz("18446744073709551616")
mul.bigz(a,a)
"340282366920938463463374607431768211456"
[edit] A native implementation
This code is more verbose than necessary, for ease of understanding.
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)
"340282366920938463463374607431768211456"
[edit] Racket
#lang racket
(define (mult A B)
(define nums
(let loop ([B B] [zeros '()])
(if (null? B)
'()
(cons (append zeros (let loop ([c 0] [A A])
(cond [(pair? A)
(define-values [q r]
(quotient/remainder
(+ c (* (car A) (car B)))
10))
(cons r (loop q (cdr A)))]
[(zero? c) '()]
[else (list c)])))
(loop (cdr B) (cons 0 zeros))))))
(let loop ([c 0] [nums nums])
(if (null? nums)
'()
(let-values ([(q r) (quotient/remainder (apply + c (map car nums)) 10)])
(cons r (loop q (filter pair? (map cdr nums))))))))
(define (number->list n)
(if (zero? n) '()
(let-values ([(q r) (quotient/remainder n 10)])
(cons r (number->list q)))))
(define 2^64 (number->list (expt 2 64)))
(for-each display (reverse (mult 2^64 2^64))) (newline)
;; for comparison
(* (expt 2 64) (expt 2 64))
;; Output:
;; 340282366920938463463374607431768211456
;; 340282366920938463463374607431768211456
[edit] REXX
[edit] version 1
This REXX version supports:
- leading signs
- decimal points
/*REXX program performs long multiplication on two numbers (without 'E')*/
numeric digits 100; d='.' /*be able to handle input numbers*/
parse arg x y . /*accept the (possible) two nums.*/
if x=='' then x=2**64 /*Not specified? Use the default*/
if y=='' then y=x /* " " " " " */
if x<0 && y<0 then sgn='-' /*only one argument is negative? */
else sgn= /*no, then the sign is positive. */
xx=x; x=strip(x,'T',d) /*remove any trailing decimal pt.*/
yy=y; y=strip(y,'T',d) /* " " " " ".*/
_=left(x,1); if _=='-' | _=='+' then x=substr(x,2) /*remove leading ±*/
_=left(y,1); if _=='-' | _=='+' then y=substr(y,2) /* " " "*/
/*[↑] above code for a Regina bug*/
/*otherwise: x=abs(x) will do it*/
dp=0; Lx=length(x); Ly=length(y) /*get the lengths of new X and Y.*/
f=pos(d,x); if f\==0 then dp= Lx-f /*calculate size of dec fraction.*/
f=pos(d,y); if f\==0 then dp=dp+Ly-f /* " " " " " */
x=space(translate(x,,d),0) /*remove decimal point, if any. */
y=space(translate(y,,d),0) /* " " " " " */
Lx=length(x); Ly=length(y) /*get the lengths of new X and Y.*/
numeric digits max(digits(),Lx+Ly)
p=0 /*the product so far. */
do j=Ly by -1 for Ly /*almost like REXX does it,but no*/
p=p+((x*substr(y,j,1))copies(0,Ly-j))
end /*j*/
say
f=length(p)-dp /*does product has enough digits?*/
if f<0 then p=copies(0,abs(f)+1)p /*Neg? Add leading 0s for INSERT*/
say ' built-in:' xx '*' yy '=' xx*yy
say 'long mult:' xx '*' yy '=' sgn ||strip(insert(d,p,length(p)-dp),'T',d)
/*stick a fork in it, we're done.*/
output when using the default inputs of: 2^64 2^64
built-in: 18446744073709551616 * 18446744073709551616 = 340282366920938463463374607431768211456 long mult: 18446744073709551616 * 18446744073709551616 = 340282366920938463463374607431768211456
output when using the inputs of: 123 -456789000
built-in: 123 * -456789000 = -56185047000 long mult: 123 * -456789000 = -56185047000
output when using the inputs of: -123.678 +456789000
built-in: -123.678 * +456789000 = -56494749942.000 long mult: -123.678 * +456789000 = -56494749942.000
[edit] version 2
/* REXX **************************************************************
* While REXX can multiply arbitrary large integers
* here is the algorithm asked for by the task description
* 13.05.2013 Walter Pachl
*********************************************************************/
cnt.=0
Numeric Digits 100
Call test 123 123
Call test 12 12
Call test 123456789012 44444444444
Call test 2**64 2**64
Call test 0 0
say cnt.0ok 'ok'
say cnt.0nok 'not ok'
Exit
test:
Parse Arg a b
soll=a*b
haben=multiply(a b)
Say 'soll =' soll
Say 'haben=' haben
If haben<>soll Then
cnt.0nok=cnt.0nok+1
Else
cnt.0ok=cnt.0ok+1
Return
multiply: Procedure
/* REXX **************************************************************
* Multiply(a b) -> a*b
*********************************************************************/
Parse Arg a b
Call s2a 'a'
Call s2a 'b'
r.=0
rim=1
r0=0
Do bi=1 To b.0
Do ai=1 To a.0
ri=ai+bi-1
p=a.ai*b.bi
Do i=ri by 1 Until p=0
s=r.i+p
r.i=s//10
p=s%10
End
rim=max(rim,i)
End
End
res=strip(a2s('r'),'L','0')
If res='' Then
res='0'
Return res
s2a:
/**********************************************************************
* copy characters of a string into a corresponding array
* digits are numbered 1 to n fron right to left
**********************************************************************/
Parse arg name
string=value(name)
lstring=length(string)
do z=1 to lstring
Call value name'.'z,substr(string,lstring-z+1,1)
End
Call value name'.0',lstring
Return
a2s:
/**********************************************************************
* turn the array of digits into a string
**********************************************************************/
call trace 'o'
Parse Arg name
ol=''
Do z=rim To 1 By -1
ol=ol||value(name'.z')
End
Return ol
Output:
soll = 15129 haben= 15129 soll = 144 haben= 144 soll = 5486968400478463649328 haben= 5486968400478463649328 soll = 340282366920938463463374607431768211456 haben= 340282366920938463463374607431768211456 soll = 0 haben= 0 5 ok 0 not ok
[edit] 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)
18446744073709551616 * 18446744073709551616 = 340282366920938463463374607431768211456 longmult(18446744073709551616, 18446744073709551616) = 340282366920938463463374607431768211456
[edit] 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.
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)))
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:
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)))
[edit] Scheme
(* (expt 2 64) (expt 2 64))
[edit] Seed7
Seed7 supports arbitrary-precision arithmetic. The library bigint.s7i defines the type bigInteger. A bigInteger is a signed integer number of unlimited size. With library support the task can be solved by using the multiplication operator *:
$ include "seed7_05.s7i";
include "bigint.s7i";
const proc: main is func
begin
writeln(2_**64 * 2_**64);
end func;
Output:
340282366920938463463374607431768211456
This task seems to prefer an inferior implementation of a long multiplication, where long numbers are stored in decimal strings. Besides type safety there are seveal other drawbacks triggered by such a representation. E.g.: In almost all cases a representation with decimal strings leads to significant lower computing speed. The multiplication example below uses the requested inferior implementation:
$ include "seed7_05.s7i";
const func string: (in string: a) * (in string: b) is func
result
var string: product is "";
local
var integer: i is 1;
var integer: j is 1;
var integer: k is 0;
var integer: carry is 0;
begin
if startsWith(a, "-") then
if startsWith(b, "-") then
product := a[2 ..] * b[2 ..];
else
product := "-" & a[2 ..] * b;
end if;
elsif startsWith(b, "-") then
product := "-" & a * b[2 ..];
else
product := "0" mult length(a) + length(b);
for i range length(a) downto 1 do
k := i + length(b);
carry := 0;
for j range length(b) downto 1 do
carry +:= (ord(a[i]) - ord('0')) * (ord(b[j]) - ord('0')) + (ord(product[k]) - ord('0'));
product @:= [k] chr(carry rem 10 + ord('0'));
carry := carry div 10;
decr(k);
end for;
product @:= [k] chr(ord(product[k]) + carry);
end for;
while startsWith(product, "0") and length(product) >= 2 do
product := product[2 ..];
end while;
end if;
end func;
const proc: main is func
begin
writeln("-18446744073709551616" * "-18446744073709551616");
end func;
The output is the same as with the superior solution.
[edit] Slate
(2 raisedTo: 64) * (2 raisedTo: 64).
[edit] Smalltalk
(2 raisedTo: 64) * (2 raisedTo: 64).
[edit] 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:
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}]
outputs
18446744073709551616 340282366920938463463374607431768211456 340282366920938463463374607431768211456
[edit] 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.
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
output:
340282366920938463463374607431768211456
[edit] Vedit macro language
This example multiplies the value on current line with the value on next line and stores result on the 3rd line.
BOL
#11 = EOL_Pos-Cur_Pos
#12 = EOL_Pos-1
Line(1)
#21 = EOL_Pos-Cur_Pos
#22 = EOL_Pos-1
EOL Ins_Newline
Ins_Char('0', COUNT, #11+#21)
#32 = Cur_Pos-1
for (#2 = 0; #2 < #21; #2++) {
Goto_Pos(#22-#2) #5 = Cur_Char - '0'
for (#1 = 0; #1 < #11; #1++) {
Goto_Pos(#12-#1) #6 = Cur_Char - '0'
#7 = #5 * #6
#3 = #1 + #2
while (#7 > 0) {
Goto_Pos(#32-#3)
#7 += Cur_Char - '0'
Ins_Char(#7%10 + '0', OVERWRITE)
#3++
#7 = #7/10
}
}
}
Sample input and output:
18446744073709551616 18446744073709551616 0340282366920938463463374607431768211456
[edit] XPL0
include c:\cxpl\stdlib;
char Two64, Product(40);
[Two64:= "18446744073709551616";
StrNMul(Two64, Two64, Product, 20);
Product(39):= Product(39)!$80; \terminate string
Text(0, Product+1); \skip leading zero
]
Output:
340282366920938463463374607431768211456
Categories:
- Programming Tasks
- Arbitrary precision
- Arithmetic operations
- Ada
- ALGOL 68
- AutoHotkey
- AWK
- BASIC
- BBC BASIC
- C
- C++
- C++ examples needing attention
- Examples needing attention
- C sharp
- F
- CoffeeScript
- Common Lisp
- D
- Dc
- Dc examples needing attention
- Euphoria
- F Sharp
- Factor
- Fortran
- Go
- Groovy
- Groovy examples needing attention
- Haskell
- Icon
- Unicon
- Icon examples needing attention
- J
- Java
- JavaScript
- Liberty BASIC
- Mathematica
- PARI/GP
- Perl
- Perl 6
- PHP
- PHP examples needing attention
- PicoLisp
- PicoLisp examples needing attention
- PL/I
- PowerShell
- Prolog
- PureBasic
- Python
- R
- Gmp
- Racket
- REXX
- Ruby
- Scala
- Scheme
- Scheme examples needing attention
- Seed7
- Slate
- Slate examples needing attention
- Smalltalk
- Smalltalk examples needing attention
- Tcl
- Ursala
- Vedit macro language
- XPL0
