Long multiplication

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

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): <lang ada>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;</lang> 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: <lang ada>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;</lang>

And finally we have the requested test application: <lang ada>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;</lang>

Output:

18446744073709551616 * 18446744073709551616 = 340282366920938463463374607431768211456

Using modular types

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

Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386

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

Works with: ALGOL 68 version Standard - no extensions to language used

Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386

<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.

AWK

Works with: gawk version 3.1.7

Works with: nawk version 20100523

Translation of: Tcl

<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

<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>

BBC BASIC

Works with: BBC BASIC for Windows

Library method: <lang bbcbasic> INSTALL @lib$+"BB4WMAPMLIB"

     MAPM_DllPath$ = @lib$+"BB4WMAPM.DLL"
     PROCMAPM_Init
     
     twoto64$ = "18446744073709551616"
     PRINT "2^64 * 2^64 = " ; FNMAPM_Multiply(twoto64$, twoto64$)</lang>

Explicit method: <lang bbcbasic> 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)</lang>

C

Doing it as if by hand. <lang c>#include <stdio.h>

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; }</lang>output<lang>340282366920938463463374607431768211456</lang>

C++

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" );

} //-------------------------------------------------------------------------------------------------- </lang>

Output:

340282366920938463463374607431768211456

Factor 1 = 9876548974569852365985574874787454878778975948
Factor 2 = 8954564845421878741168741154541897945138974567
Product: = 88440198241770705041777453160463400993104404280916080859287340887463980926235972531076714516

Version 2

This example is incorrect. Please fix the code and remove this message.

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#

This example is incorrect. Please fix the code and remove this message.

System.Numerics.BigInteger was added with C# 4.

Works with: C# version 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

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 </lang>

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>void main() {

   import std.stdio, std.bigint;

   writeln(2.BigInt ^^ 64 * 2.BigInt ^^ 64);

}</lang>

Output:

340282366920938463463374607431768211456

Long multiplication, same output:

Translation of: JavaScript

<lang d>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;

}</lang>

Dc

This example is incorrect. Please fix the code and remove this message.

Since Dc has arbitrary precision built-in, the task is no different than a normal multiplication: <lang Dc>2 64^ 2 64^ *p</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#

This example is incorrect. Please fix the code and remove this message.

<lang F#>> let X = 2I ** 64 * 2I ** 64 ;;

val X : System.Numerics.BigInteger = 340282366920938463463374607431768211456 </lang>

Factor

<lang 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 ;</lang>

<lang factor>( scratchpad ) 2 64 ^ dup longmult . 340282366920938463463374607431768211456 ( scratchpad ) 2 64 ^ dup * . 340282366920938463463374607431768211456</lang>

Fortran

Works with: Fortran version 95 and later

<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>

Go

<lang 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))

}</lang> Output:

340282366920938463463374607431768211456

Groovy

This example is incorrect. Please fix the code and remove this message.

<lang groovy> println 2**64 * 2**64 </lang>

340282366920938463463374607431768211456

Haskell

<lang 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)</lang>

Output:

340282366920938463463374607431768211456

Icon and Unicon

This example is incorrect. Please fix the code and remove this message.

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:

digits translates a number to a corresponding list of digits;

polymult (multiplies polynomials): ref. [1]

buildDecimal (translates a list of decimal digits to the corresponding extended precision number):

Java

This example is incorrect. Please fix the code and remove this message.

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> '[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 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$

   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

</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 $3^{(10^{7})}\times 3^{(10^{7})}$ in about a second; the final result has 9542426 digits; result omitted for obvious reasons.

PARI/GP

<lang parigp>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")</lang> Output:

%1 = "340282366920938463463374607431768211456"

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 two implementations (Pugs and Niecza) support already:

niecza> 18446744073709551616 * 18446744073709551616
340282366920938463463374607431768211456

PHP

This example is incorrect. Please fix the code and remove this message.

<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

This example is incorrect. Please fix the code and remove this message.

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

This example is incorrect. Please fix the code and remove this message.

PureBasic

This example is incorrect. Please fix the code and remove this message.

Works with: PureBasic version 4.41

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>

Works with: Python version 3.0

Translation of: Perl

<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:

Translation of: Haskell

<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

Library: 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"

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

</lang>

REXX

version 1

This REXX version supports:

leading signs
decimal points

<lang rexx>/*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.*/</lang>

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

version 2

<lang rexx>/* REXX **************************************************************

While REXX can multiply arbitrary lage 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</lang>

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

Ruby

Translation of: Tcl

<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

Works with: Scala version 2.8

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

This example is incorrect. Please fix the code and remove this message.

Seed7

This example is incorrect. Please fix the code and remove this message.

<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

This example is incorrect. Please fix the code and remove this message.

<lang slate>(2 raisedTo: 64) * (2 raisedTo: 64).</lang>

Smalltalk

This example is incorrect. Please fix the code and remove this message.

<lang smalltalk>(2 raisedTo: 64) * (2 raisedTo: 64).</lang>

Tcl

Works with: Tcl version 8.5

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

Vedit macro language

This example multiplies the value on current line with the value on next line and stores result on the 3rd line. <lang vedit>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

       }
   }

} </lang> Sample input and output:

18446744073709551616
18446744073709551616
0340282366920938463463374607431768211456

XPL0

<lang 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 ]</lang>

Output:

340282366920938463463374607431768211456