Long multiplication: Difference between revisions

From Rosetta Code
Content added Content deleted
(add AWK)
(→‎{{header|Python}}: (Note that Python comes with arbitrary length integers).)
Line 603: Line 603:


=={{header|Python}}==
=={{header|Python}}==
(Note that Python comes with arbitrary length integers).

{{works with|Python|3.0}}
{{works with|Python|3.0}}
{{trans|Perl}}
{{trans|Perl}}

Revision as of 17:39, 4 June 2009

Task
Long multiplication
You are encouraged to solve this task according to the task description, using any language you may know.

In this task, explicitly implement long multiplication. This is one possible approach to arbitrary-precision integer algebra.

For output, display the result of 2^64 * 2^64. The decimal representation of 2^64 is: 

18446744073709551616

The output of 2^64 * 2^64 is 2^128, and that is: 

340282366920938463463374607431768211456

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 algol>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 algol>MODE DIGIT = INT; MODE INTEGER = FLEX[0]DIGIT; # an arbitary number of digits #

  1. "digits" are stored in digit base ten, but 10000 & 2**n (inc hex) can be used #

INT digit base = 1000;

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

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

);

  1. initialise/cast a INTEGER from an LONG INT #

OP INTEGERINIT = (LONG INT number)INTEGER: INTEGERINIT LENG number;

  1. initialise/cast a INTEGER from an INT #

OP INTEGERINIT = (INT number)INTEGER: INTEGERINIT LENG LENG number;

  1. 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);
  1. 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;
  1. 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 | "-" | "+" );
  1. finally return the semi-normalised result #
 IF leading zeros = UPB out THEN "0" ELSE sign + out[leading zeros+1:] FI

);</lang> <lang algol>################################################################

  1. Finally Define the required INTEGER multiplication OPerator. #

OP * = (INTEGER a, b)INTEGER:(

  1. 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;
  1. 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 algol># 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.0
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 = 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) { 

   split(x, a_x, //)
   return reverse_array(a_x, a)

}

function reverse_array(a,b) { 

   len = length_array(a)
   for (i in a) {
       b[1+len-i] = a[i]
   }
   return len

}

  1. work around a bug in gawk 3.1.6 where length(array) doesn't work
  2. if the array is a function parameter

function length_array(a) { 

   len = asort(a, x)
   return len

}

function append_array(a, value) { 

   len = length_array(a)
   a[++len] = value

}

function reverse_join(a) { 

   len = length_array(a)
   str = ""
   for (i=len; i>=1; i--) {
       str = str a[i]
   }
   return str

}

function print_array(a) { 

   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

C

<lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <string.h>
  1. define X(V,L,I) ( ((I)<(L)) ? (V[(L)-(I)-1]-'0') : 0)
  2. define MIN(A,B) ( ((A)<(B)) ? (A) : (B) )

char *longmolt(char *a, char *b) { 

 int n, m, T;
 char *r = NULL;
 int i, j, k, p;
 int *C, *R;
 
 n = strlen(a);
 m = strlen(b);
 T = n+m+2;

 r = malloc(T+1);
 R = malloc(T*sizeof(int));
 C = malloc((n+1)*(m+1)*sizeof(int));

 memset(r, '0', T);
 memset(C, 0, (n+1)*(m+1)*sizeof(int));
 r[T] = 0;

 for(i=0; i<(m+1); i++)
 {
   C[i*(n+1)] = (X(b,m,i) * X(a,n,0));
   for(j=1; j<(n+1); j++)
   {
     C[i*(n+1)+j] = (X(b,m,i) * X(a,n,j) + C[i*(n+1)+j-1] / 10);
   }
 }
 for(k=0; k < T; k++)
 {
   R[k] = 0;
   for(j=0; j < MIN(k,m+1) ; j++)
   {
     i = k-j-1;
     if ( (i>n) || (i<0) ) { continue; }
     R[k] += (C[j*(n+1)+i]%10);
   }
   R[k] += ((k-1)<0 ) ? 0 : (R[k-1]/10);
 }
 for(k=T; k>0; k--) r[k-1] = R[T-k+1]%10 + '0';
 free(C); free(R);
 return r;

}</lang>

<lang c>/* using */ char *n1 = "18446744073709551616"; int main() { 

 char *res;
 int lz;

 /* printf("%s * %s = 340282366920938463463374607431768211456\n", n1, n1); */
 res = longmolt(n1,n1);
 for(lz=0; (lz < strlen(res)) && (res[lz]=='0') ;lz++) ; 
 printf("%s * %s = %s\n", n1, n1, res+lz);
 free(res);
 return 0;

}</lang>

The code does not handle negative intergers, nor there are proper error checks. It is more or less the implementation of the way a human being multiplies two integers. The numbers are stored as strings.

Using GMP (GNU Multi Precision library)

Library: GMP

<lang c>#include <stdio.h>

  1. include <gmp.h>

int main() { 

 mpz_t z1, z2, zr;

 mpz_init(z2);   mpz_init(zr);
 mpz_init_set_str(z1, "18446744073709551616", 10);
 mpz_set(z2, z1);
 mpz_mul(zr, z1, z2);
 mpz_out_str(stdout, 10, zr);
 printf("\n");
 return 0;

}</lang>

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>

Haskell

<lang haskell>digits :: Integer -> [Integer] digits = map (fromIntegral.digitToInt) . show

lZZ = inits $ repeat 0

table f = map . flip (map . f)

polymul = ((map sum . transpose . zipWith (++) lZZ) .) . table (*)

longmult = (foldl1 ((+) . (10 *)) .) . (. digits) . polymul . digits</lang> Output: <lang haskell>*Main> (2^64) `longmult` (2^64) 340282366920938463463374607431768211456</lang>

J

   (([+10x*])/@|.@(,.&.":@[+//.@(*/),.&.":@]))/ ,~2x^64
340282366920938463463374607431768211456
  • digits: ,.&.": y
   ,.&.": 123
1 2 3
  • polynomial multiplication: x (+//.@(*/)) y ref. [1]
   1 2 3 (+//.@(*/)) 1 2 3
1 4 10 12 9
  • building the decimal result: ([+10x*])/|. y
   ([+10x*])/|. 1 4 10 12 9
15129

or using the primitive dyad #. instead of ([+10x*])/@|. 

   (10x #.,.&.":@[+//.@(*/),.&.":@])/ ,~2x^64
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>


Perl

<lang perl>#!/usr/bin/perl -w use strict;

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

Python

(Note that Python comes with arbitrary length integers).

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>

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
Retrieved from "https://rosettacode.org/wiki/Long_multiplication?oldid=101190"