Long multiplication

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

Explicitly implement   long multiplication.

This is one possible approach to arbitrary-precision integer algebra.


For output, display the result of   264 * 264.

Optionally, verify your result against builtin arbitrary precision support.

The decimal representation of   264   is:

18,446,744,073,709,551,616

The output of   264 * 264   is   2128,   and is:

340,282,366,920,938,463,463,374,607,431,768,211,456



11l

Translation of: Python
F add_with_carry(&result, =addend, =addendpos)
   L
      L result.len < addendpos + 1
         result.append(‘0’)
      V addend_result = String(Int(addend) + Int(result[addendpos]))
      V addend_digits = Array(addend_result)
      result[addendpos] = addend_digits.pop()

      I addend_digits.empty
         L.break

      addend = addend_digits.pop()
      addendpos++

F longhand_multiplication(multiplicand, multiplier)
   [Char] result
   L(multiplicand_digit) reversed(multiplicand)
      V multiplicand_offset = L.index
      L(multiplier_digit) reversed(multiplier)
         V multiplier_offset = L.index + multiplicand_offset
         V multiplication_result = String(Int(multiplicand_digit) * Int(multiplier_digit))

         L(result_digit_addend) reversed(multiplication_result)
            V addend_offset = L.index + multiplier_offset
            add_with_carry(&result, result_digit_addend, addend_offset)

   result.reverse()
   R result.join(‘’)

V sixtyfour = ‘18446744073709551616’
print(longhand_multiplication(sixtyfour, sixtyfour))
Output:
340282366920938463463374607431768211456

360 Assembly

For maximum compatibility, we use only the basic 370 instruction set (use of MVCL). Pseudo-macro instruction XPRNT can be replaced by a WTO.

LONGINT  CSECT
         USING  LONGINT,R13
SAVEAREA B      PROLOG-SAVEAREA(R15)
         DC     17F'0'
         DC     CL8'LONGINT'
PROLOG   STM    R14,R12,12(R13)
         ST     R13,4(R15)
         ST     R15,8(R13)
         LR     R13,R15
         MVC    XX(1),=C'1'
         MVC    LENXX,=H'1'        xx=1
         LA     R2,64
LOOPII   ST     R2,RLOOPII         do for 64
         MVC    X-2(LL+2),XX-2     x=xx
         MVC    Y(1),=C'2'
         MVC    LENY,=H'1'         y=2
         BAL    R14,LONGMULT
         MVC    XX-2(LL+2),Z-2     xx=longmult(xx,2)   xx=xx*2
         L      R2,RLOOPII
ELOOPII  BCT    R2,LOOPII          loop
         MVC    X-2(LL+2),XX-2
         MVC    Y-2(LL+2),XX-2
         BAL    R14,LONGMULT
         MVC    YY-2(LL+2),Z-2     yy=longmult(xx,xx)  yy=xx*xx
         XPRNT  XX,LL              output xx
         XPRNT  YY,LL              output yy
RETURN   L      R13,4(0,R13)       epilog
         LM     R14,R12,12(R13)
         XR     R15,R15            set return code
         BR     R14                return to caller
RLOOPII  DS     F
*        
LONGMULT EQU    *                  function longmult z=(x,y)
         MVC    LENSHIFT,=H'0'     shift=''
         MVC    LENZ,=H'0'         z=''
         LH     R6,LENX     
         LA     R6,1(R6)           from lenx     
         XR     R8,R8
         BCTR   R8,0               by -1
         LA     R9,0               to 1
LOOPI    BXLE   R6,R8,ELOOPI       do i=lenx to 1 by -1
         LA     R2,X
         AR     R2,R6              +i
         BCTR   R2,0
         MVC    CI,0(R2)           ci=substr(x,i,1)
         IC     R0,CI              ni=integer(ci)
         N      R0,=X'0000000F'
         STH    R0,NI
         MVC    LENT,=H'0'         t=''
         SR     R0,R0
         STH    R0,CARRY           carry=0
         LH     R7,LENY
         LA     R7,1(R7)           from lenx     
         XR     R10,R10
         BCTR   R10,0              by -1
         LA     R11,0              to 1
LOOPJ1   BXLE   R7,R10,ELOOPJ1     do j=leny to 1 by -1 
         LA     R2,Y
         AR     R2,R7              +j
         BCTR   R2,0
         MVC    CJ,0(R2)           cj=substr(y,j,1)
         IC     R0,CJ
         N      R0,=X'0000000F'
         STH    R0,NJ              nj=integer(cj)
         LH     R2,NI
         MH     R2,NJ
         AH     R2,CARRY
         STH    R2,NKR             nkr=ni*nj+carry
         LH     R2,NKR
         LA     R1,10
         SRDA   R2,32
         DR     R2,R1
         STH    R2,NK              nk=nkr//10
         STH    R3,CARRY           carry=nkr/10
         LH     R2,NK
         O      R2,=X'000000F0'
         STC    R2,CK              ck=string(nk)
         MVC    TEMP,T
         MVC    T(1),CK
         MVC    T+1(LL-1),TEMP
         LH     R2,LENT
         LA     R2,1(R2)
         STH    R2,LENT            t=ck!!t
         B      LOOPJ1             next j
ELOOPJ1  EQU    *
         LH     R2,CARRY
         O      R2,=X'000000F0'
         STC    R2,CK              ck=string(carry)
         MVC    TEMP,T
         MVC    T(1),CK
         MVC    T+1(LL-1),TEMP
         LH     R2,LENT
         LA     R2,1(R2)
         STH    R2,LENT            t=ck!!t
         LA     R2,T
         AH     R2,LENT
         LH     R3,LENSHIFT
         LA     R4,SHIFT
         LH     R5,LENSHIFT
         MVCL   R2,R4
         LH     R2,LENT
         AH     R2,LENSHIFT
         STH    R2,LENT            t=t!!shift
IF1      LH     R4,LENZ
         CH     R4,LENT            if lenz>lent
         BNH    ELSE1
         LH     R2,LENZ            then
         LA     R2,1(R2)
         STH    R2,L               l=lenz+1
         B      EIF1
ELSE1    LH     R2,LENT            else 
         LA     R2,1(R2)
         STH    R2,L               l=lent+1
EIF1     EQU    *
         MVI    TEMP,C'0'          to
         MVC    TEMP+1(LL-1),TEMP
         LA     R2,TEMP
         AH     R2,L
         SH     R2,LENZ
         LH     R3,LENZ
         LA     R4,Z
         LH     R5,LENZ
         MVCL   R2,R4
         MVC    LENZ,L
         MVC    Z,TEMP             z=right(z,l,'0')
         MVI    TEMP,C'0'          to
         MVC    TEMP+1(LL-1),TEMP
         LA     R2,TEMP
         AH     R2,L
         SH     R2,LENT
         LH     R3,LENT
         LA     R4,T
         LH     R5,LENT
         MVCL   R2,R4
         MVC    LENT,L
         MVC    T,TEMP             t=right(t,l,'0')
         MVC    LENW,=H'0'         w=''
         SR     R0,R0
         STH    R0,CARRY           carry=0
         LH     R7,L
         LA     R7,1(R7)           from l
         XR     R10,R10
         BCTR   R10,0              by -1
         LA     R11,0              to 1
LOOPJ2   BXLE   R7,R10,ELOOPJ2     do j=l to 1 by -1
         LA     R2,Z
         AR     R2,R7              +j
         BCTR   R2,0
         MVC    CZ,0(R2)           cz=substr(z,j,1)
         IC     R0,CZ
         N      R0,=X'0000000F'
         STH    R0,NZ              nz=integer(cz)
         LA     R2,T
         AR     R2,R7              -j
         BCTR   R2,0
         MVC    CT,0(R2)           ct=substr(t,j,1)
         IC     R0,CT
         N      R0,=X'0000000F'
         STH    R0,NT              nt=integer(ct)
         LH     R2,NZ
         AH     R2,NT
         AH     R2,CARRY
         STH    R2,NKR             nkr=nz+nt+carry
         LH     R2,NKR
         LA     R1,10
         SRDA   R2,32
         DR     R2,R1
         STH    R2,NK
         STH    R3,CARRY           nk=nkr//10; carry=nkr/10
         LH     R2,NK
         O      R2,=X'000000F0'
         STC    R2,CK              ck=string(nk)
         MVC    TEMP,W
         MVC    W(1),CK
         MVC    W+1(LL-1),TEMP
         LH     R2,LENW
         LA     R2,1(R2)
         STH    R2,LENW            w=ck!!w
         B      LOOPJ2             next j
ELOOPJ2  EQU    *
         LH     R2,CARRY
         O      R2,=X'000000F0'
         STC    R2,CK             ck=string(carry)
         MVC    Z(1),CK
         MVC    Z+1(LL-1),W
         LH     R2,LENW
         LA     R2,1(R2)
         STH    R2,LENZ            z=ck!!w
         LA     R7,0               from 1
         LA     R10,1              by 1
         LH     R11,LENZ           to lenz
LOOPJ3   BXH    R7,R10,ELOOPJ3     do j=1 to lenz
         LA     R2,Z
         AR     R2,R7              j
         BCTR   R2,0
         MVC    ZJ(1),0(R2)        zj=substr(z,j,1)
         CLI    ZJ,C'0'            if zj^='0'
         BNE    ELOOPJ3            then leave j
         B      LOOPJ3             next j
ELOOPJ3  EQU    *
IF2      CH     R7,LENZ            if j>lenz
         BNH    EIF2
         LH     R7,LENZ            then j=lenz
EIF2     EQU    *
         LA     R2,TEMP            to
         LH     R3,LENZ
         SR     R3,R7              -j
         LA     R3,1(R3)
         STH    R3,LENTEMP
         LA     R4,Z               from
         AR     R4,R7              +j
         BCTR   R4,0
         LR     R5,R3
         MVCL   R2,R4
         MVC    Z-2(LL+2),TEMP-2   z=substr(z,j)
         LA     R2,SHIFT
         AH     R2,LENSHIFT
         MVI    0(R2),C'0'
         LH     R3,LENSHIFT
         LA     R3,1(R3)
         STH    R3,LENSHIFT        shift=shift!!'0'
         MVC    TEMP,Z
         LA     R2,TEMP
         AH     R2,LENZ
         MVC    0(2,R2),=C'  '
         B      LOOPI              next i
ELOOPI   EQU    *
         MVI    TEMP,C' '
         LA     R2,Z
         AH     R2,LENZ
         LH     R3,=AL2(LL)
         SH     R3,LENZ
         LA     R4,TEMP
         LH     R5,=H'1'
         ICM    R5,8,=C' '
         MVCL   R2,R4              z=clean(z)
         BR     R14                end function longmult
*
L        DS     H
NI       DS     H
NJ       DS     H
NK       DS     H
NZ       DS     H
NT       DS     H
CARRY    DS     H
NKR      DS     H
CI       DS     CL1
CJ       DS     CL1
CZ       DS     CL1
CT       DS     CL1
CK       DS     CL1
ZJ       DS     CL1
LENXX    DS     H
XX       DS     CL94
LENYY    DS     H
YY       DS     CL94
LENX     DS     H
X        DS     CL94
LENY     DS     H
Y        DS     CL94
LENZ     DS     H
Z        DS     CL94
LENT     DS     H
T        DS     CL94
LENW     DS     H
W        DS     CL94
LENSHIFT DS     H
SHIFT    DS     CL94
LENTEMP  DS     H
TEMP     DS     CL94
LL       EQU    94
         YREGS  
         END    LONGINT
Output:
18446744073709551616
340282366920938463463374607431768211456

AArch64 Assembly

Works with: as version Raspberry Pi 3B version Buster 64 bits
or android 64 bits with application Termux
/* ARM assembly AARCH64 Raspberry PI 3B */
/*  program longmulti64.s   */
/* REMARK : this program use factors unsigned to 2 power 127
              and the result is less than 2 power 255 */
/************************************/
/* Constantes                       */
/************************************/
/* for this file see task include a file in language AArch64 assembly*/
.include "../includeConstantesARM64.inc" 
.equ BUFFERSIZE,   100

/***********************************************/
/* structures                                  */
/**********************************************/
/* Définition multi128 */
    .struct  0
multi128_N1:                      //   63-0
    .struct  multi128_N1 + 8  
multi128_N2:                      //  127-64
    .struct  multi128_N2 + 8 
multi128_N3:                      //  128-191
    .struct  multi128_N3 + 8 
multi128_N4:                      //  192-255
    .struct  multi128_N4 + 8 
multi128_end:
/*********************************/
/* Initialized data              */
/*********************************/
.data
szMessFactor:         .asciz "Factor = "
szMessResult:         .asciz "Result = "
szMessStart:          .asciz "Program 64 bits start.\n"
szCarriageReturn:     .asciz "\n"

i128test1:            .quad  0,1,0,0   // 2 power 64
/*********************************/
/* UnInitialized data            */
/*********************************/
.bss
sZoneConv:             .skip BUFFERSIZE   // conversion buffer
i128Result1:           .skip multi128_end 

/*********************************/
/*  code section                 */
/*********************************/
.text
.global main 
main:                            // entry of program 
    ldr x0,qAdrszMessStart
    bl affichageMess
    ldr x0,qAdri128test1         // origin number
    ldr x1,qAdrsZoneConv
    mov x2,#BUFFERSIZE
    bl convertMultiForString     // convert multi number to string
    mov x2,x0                    // insert conversion in message
    mov x0,#3                    // string number to display
    ldr x1,qAdrszMessFactor
    ldr x3,qAdrszCarriageReturn
    bl displayStrings            // display message
 
    // multiplication
    ldr x0,qAdri128test1         // factor 1
    ldr x1,qAdri128test1         // factor 2
    ldr x2,qAdri128Result1       // result 
    bl multiplierMulti128
    ldr x0,qAdri128Result1
    ldr x1,qAdrsZoneConv
    mov x2,#BUFFERSIZE
    bl convertMultiForString     // conversion multi to string
    mov x2,x0                    // insert conversion in message
    mov x0,#3                    // number string to display
    ldr x1,qAdrszMessResult
    ldr x3,qAdrszCarriageReturn
    bl displayStrings            // display message
 
100:                             // standard end of the program 
    mov x0, #0                   // return code
    mov x8,EXIT 
    svc #0                       // perform the system call
    
qAdrszCarriageReturn:        .quad szCarriageReturn
qAdrsZoneConv:               .quad sZoneConv
qAdri128test1:               .quad i128test1
qAdri128Result1:             .quad i128Result1
qAdrszMessResult:            .quad szMessResult
qAdrszMessFactor:            .quad szMessFactor
qAdrszMessStart:             .quad szMessStart
/***************************************************/
/*  multiplication multi128 by multi128             */
/***************************************************/
// x0 contains address multi128 1
// x1 contains address multi128 2
// x2 contains address result multi128
// x0 return address result (= x2)
multiplierMulti128:
    stp x1,lr,[sp,-16]!          // save  registers 
    mov x9,x0                    // factor 1
    mov x10,x1                   // factor 2
    mov x7,x2                    // address result
    mov x6,#3                    // multi128 size 
1:
    str xzr,[x7,x6,lsl #3]       // init result
    subs x6,x6,#1
    bge 1b
    mov x5,#0                    // indice loop 1
2:                               // loop items factor 1
    ldr x0,[x9,x5,lsl #3]        // load a item
    mov x4,#0
    mov x8,#0
3:                               // loop item factor 2
    add x6,x4,x5                 // compute result indice
 
    ldr x1,[x10,x4,lsl #3]       // load a item factor 2
    mul x2,x1,x0                 // multiply low 64 bits 
    umulh x3,x1,x0               // multiply high 64 bits
    ldr x1,[x7,x6,lsl #3]        // load previous item of result
    adds x1,x1,x2                // add low part result multiplication
    mov x11,1
    csel x2,x11,xzr,cs
    adds  x1,x1,x8               // add high part precedente
    adc  x8,x3,x2                // new high part with retenue
    str x1,[x7,x6,lsl #3]        // store the sum in result
  
    add x4,x4,#1
    cmp x4,#3
    blt 3b                       // and loop  2
    cmp x8,#0                    // high part ?
    beq 5f
    add x6,x6,#1
    cmp x6,#2                    // on last item ?
    ble 4f
    adr x0,szMessErrOverflow     // yes -> overflow
    bl affichageMess
    mov x0,#0                    // return 0
    b 100f
4:
    str x8,[x7,x6,lsl #3]        // no store high part in next item
5:
    add x5,x5,#1 
    cmp x5,#3
    blt 2b                       // and loop 1
    mov x0,x7

100:
    ldp x1,lr,[sp],16            // restaur  registers
    ret 
szMessErrOverflow:        .asciz "\033[31mOverflow !!\033[0m \n"
.align 4

/***************************************************/
/*   conversion multi128 unsigned to string    */
/***************************************************/
// x0 contains address multi128
// x1 contains address buffer
// x2 contains buffer length
convertMultiForString:
    stp x1,lr,[sp,-16]!          // save  registers 
    stp x2,x3,[sp,-16]!          // save  registers 
    stp x4,x5,[sp,-16]!          // save  registers 
    sub sp,sp,#multi128_end      // reserve place to stack
    mov fp,sp                    // init address to quotient
    mov x5,x1                    // save address buffer
    mov x3,#0                    // init indice
1:
    ldr x4,[x0,x3,lsl #3]        // load one part of number
    str x4,[fp,x3,lsl #3]        // copy part on stack
    add x3,x3,#1
    cmp x3,#4
    blt 1b
 
2:
    strb wzr,[x5,x2]             // store final 0 in buffer
    sub x4,x2,#1                 // end number storage
3:
    mov x0,fp
    mov x1,#10
    bl calculerModuloMultiEntier // compute modulo 10
    add x0,x0,#0x30              // convert result to character
    strb w0,[x5,x4]              // store character on buffer
    subs x4,x4,#1                // 
    blt 99f                      //  buffer too low
    ldr x0,[fp,#multi128_N1]     // test if quotient = zero
    cmp x0,#0
    bne 3b
    ldr x0,[fp,#multi128_N2]
    cmp x0,#0
    bne 3b
    ldr x0,[fp,#multi128_N3]
    cmp x0,#0
    bne 3b
    ldr x0,[fp,#multi128_N4]
    cmp x0,#0
    bne 3b
 

    add x0,x5,x4                 // return begin number in buffer
    add x0,x0,#1            
    b 100f
99:                              // display error if buffer est toop low
    adr x0,szMessErrBuffer
    bl affichageMess
    mov x0,#-1
100:
    add sp,sp,#multi128_end      // stack alignement
    ldp x4,x5,[sp],16        // restaur  registers 
    ldp x2,x3,[sp],16        // restaur  registers 
    ldp x1,lr,[sp],16        // restaur  registers
    ret 
szMessErrBuffer:        .asciz "\033[31mBuffer de conversion trop petit !!\033[0m \n"
.align 4
/***************************************************/
/*    modulo  compute   unsigned                   */
/***************************************************/
// x0 contains address multi128
// x1 contains modulo (positive)
// x0 return  modulo
// ATTENTION : le multientier origine est modifié et contient le quotient
calculerModuloMultiEntier: // INFO: calculerModuloMultiEntier
    stp x1,lr,[sp,-16]!    // save  registers 
    stp x2,x3,[sp,-16]!    // save  registers 
    stp x4,x5,[sp,-16]!    // save  registers 
    cmp x1,#0
    ble 99f
    mov x4,x1              // save modulo
    mov x3,#3
    mov x5,x0              // multi128 address
    ldr x0,[x5,x3,lsl 3]   // load last part of number in low part of 128 bits
    mov x1,#0              // init higt part 128 bits
1:
    cmp x3,#0              // end part ?
    ble 2f
    mov x2,x4              // modulo
    bl division64R         // divide x0,x1 by x2 in x0,x1 and remainder in x2
    str x0,[x5,x3,lsl #3]  // store result part low  
    sub x3,x3,#1           // other part ?
    ldr x0,[x5,x3,lsl #3]  // load prev part
    mov x1,x2              // store remainder on high part of 128 bits           
    b 1b
2:
    mov x2,x4              // modulo
    bl division64R
    str x0,[x5]            // stockage dans le 1er chunk
    mov x0,x2              // return remainder
    b 100f
99:
    adr x0,szMessNegatif
    bl affichageMess
    mov x0,#-1
100:                       // fin standard de la fonction  
    ldp x4,x5,[sp],16      // restaur  registers 
    ldp x2,x3,[sp],16      // restaur  registers 
    ldp x1,lr,[sp],16      // restaur  registers
    ret 
szMessNegatif:      .asciz "\033[31mLe diviseur doit être positif !\033[0m\n"
.align 4
/***************************************************/
/*   division 128 bits number in 2 registers by 64 bits number */
/***************************************************/
/* x0 contains  dividende  low part */
/* x1 contains  dividende  high part */
/* x2 contains divisor             */
/* x0 return quotient  low part  */
/* x1 return quotient  high part  */
/* x2 return remainder             */
division64R:
    stp x3,lr,[sp,-16]!  // save  registers 
    stp x4,x5,[sp,-16]!  // save  registers 
    stp x6,x7,[sp,-16]!  // save  registers 
    stp x8,x9,[sp,-16]!  // save  registers 
    mov x6,#0            // init high high part of remainder !!
                         // x1 = high part of number in high part of remainder
    mov x7,x0            // low part of number in low part of remainder
    mov x3,#0            // init high part quotient
    mov x4,#0            // init low part quotient
    mov x5,#64
1:                       // begin loop
    lsl x6,x6,#1         // left shift high high part of remainder
    cmp x1,0             // if negative ie bit 63 = 1
    orr x8,x6,1
    csel x6,x8,x6,lt     // add left bit  high part on high high part
    lsl x1,x1,#1         // left shift high part of remainder 
    cmp x7,0
    orr x8,x1,1
    csel x1,x8,x1,lt     // add left bit  low part on high part
    lsl x7,x7,#1         // left shift low part of remainder 
    cmp x4,0
    lsl x4,x4,#1         // left shift low part quotient 
    lsl x3,x3,#1         // left shift high part quotient
    orr x8,x3,1
    csel x3,x8,x3,lt     // add left bit low part on high part
                         // sub divisor to high part remainder
    subs x1,x1,x2
    sbcs  x6,x6,xzr      // sub restr.quad (retenue in french)
    bmi 2f               // result negative ?
                         // positive or equal
    orr x4,x4,#1         // right bit quotient  to 1
    b 3f
2:                       // negative
    orr x4,x4,xzr        // right bit quotient to 0 
    adds x1,x1,x2        // and restaure the remainder to precedent value
    adc  x6,x6,xzr       // and restr.quad
3:
    subs x5,x5,#1        // decrement indice 
    bgt 1b               // and loop
    mov x0,x4            // low part quotient 
    mov x2,x1            // remainder
    mov x1,x3            // high part quotient
100:
    ldp x8,x9,[sp],16    // restaur  registers 
    ldp x6,x7,[sp],16    // restaur  registers 
    ldp x4,x5,[sp],16    // restaur  registers 
    ldp x3,lr,[sp],16    // restaur  registers
    ret 
/***************************************************/
/*   display multi strings                         */
/*   new version 24/05/2023                        */
/***************************************************/
/* x0  contains number strings address */
/* x1 address string1 */
/* x2 address string2 */
/* x3 address string3 */
/* x4 address string4 */
/* x5 address string5 */
/* x6 address string5 */
displayStrings:            // INFO:  displayStrings
    stp x7,lr,[sp,-16]!          // save  registers 
    stp x2,fp,[sp,-16]!          // save  registers 
    add fp,sp,#32          // save paraméters address (4 registers saved * 8 bytes)
    mov x7,x0              // save strings number
    cmp x7,#0              // 0 string -> end
    ble 100f
    mov x0,x1              // string 1
    bl affichageMess
    cmp x7,#1              // number > 1
    ble 100f
    mov x0,x2
    bl affichageMess
    cmp x7,#2
    ble 100f
    mov x0,x3
    bl affichageMess
    cmp x7,#3
    ble 100f
    mov x0,x4
    bl affichageMess
    cmp x7,#4
    ble 100f
    mov x0,x5
    bl affichageMess
    cmp x7,#5
    ble 100f
    mov x0,x6
    bl affichageMess
100:
    ldp x2,fp,[sp],16        // restaur  registers 
    ldp x7,lr,[sp],16            // restaur  registers
    ret
/***************************************************/
/*      ROUTINES INCLUDE                           */
/***************************************************/
/* for this file see task include a file in language AArch64 assembly*/
.include "../includeARM64.inc"
Output:
Program 64 bits start.
Factor = 18446744073709551616
Result = 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):

package Long_Multiplication is
   type Number (<>) is private;

   Zero : constant Number;
   One  : constant Number;

   function Value (Item : in String) return Number;
   function Image (Item : in Number) return String;

   overriding
   function "=" (Left, Right : in Number) return Boolean;

   function "+" (Left, Right : in Number) return Number;
   function "*" (Left, Right : in Number) return Number;

   function Trim (Item : in Number) return Number;
private
   Bits : constant := 16;
   Base : constant := 2 ** Bits;

   type Accumulated_Value is range 0 .. (Base - 1) * Base;
   subtype Digit is Accumulated_Value range 0 .. Base - 1;

   type Number is array (Natural range <>) of Digit;
   for Number'Component_Size use Bits; -- or pragma Pack (Number);

   Zero : constant Number := (1 .. 0 => 0);
   One  : constant Number := (0 => 1);

   procedure Divide (Dividend  : in     Number;
                     Divisor   : in     Digit;
                     Result    :    out Number;
                     Remainder :    out Digit);
end Long_Multiplication;

Some of the operations declared above are useful helper operations for the conversion of numbers to and from base 10 digit strings.

Then we implement the operations:

package body Long_Multiplication is
   function Value (Item : in String) return Number is
      subtype Base_Ten_Digit is Digit range 0 .. 9;
      Ten : constant Number := (0 => 10);
   begin
      case Item'Length is
         when 0 =>
            raise Constraint_Error;
         when 1 =>
            return (0 => Base_Ten_Digit'Value (Item));
         when others =>
            return (0 => Base_Ten_Digit'Value (Item (Item'Last .. Item'Last)))
              + Ten * Value (Item (Item'First .. Item'Last - 1));
      end case;
   end Value;

   function Image (Item : in Number) return String is
      Base_Ten  : constant array (Digit range 0 .. 9) of String (1 .. 1) :=
                    ("0", "1", "2", "3", "4", "5", "6", "7", "8", "9");
      Result    : Number (0 .. Item'Last);
      Remainder : Digit;
   begin
      if Item = Zero then
         return "0";
      else
         Divide (Dividend  => Item,
                 Divisor   => 10,
                 Result    => Result,
                 Remainder => Remainder);

         if Result = Zero then
            return Base_Ten (Remainder);
         else
            return Image (Trim (Result)) & Base_Ten (Remainder);
         end if;
      end if;
   end Image;

   overriding
   function "=" (Left, Right : in Number) return Boolean is
   begin
      for Position in Integer'Min (Left'First, Right'First) ..
                      Integer'Max (Left'Last,  Right'Last) loop
         if Position in Left'Range and Position in Right'Range then
            if Left (Position) /= Right (Position) then
               return False;
            end if;
         elsif Position in Left'Range then
            if Left (Position) /= 0 then
               return False;
            end if;
         elsif Position in Right'Range then
            if Right (Position) /= 0 then
               return False;
            end if;
         else
            raise Program_Error;
         end if;
      end loop;

      return True;
   end "=";

   function "+" (Left, Right : in Number) return Number is
      Result      : Number (Integer'Min (Left'First, Right'First) ..
                            Integer'Max (Left'Last , Right'Last) + 1);
      Accumulator : Accumulated_Value := 0;
      Used        : Integer := Integer'First;
   begin
      for Position in Result'Range loop
         if Position in Left'Range then
            Accumulator := Accumulator + Left (Position);
         end if;

         if Position in Right'Range then
            Accumulator := Accumulator + Right (Position);
         end if;

         Result (Position) := Accumulator mod Base;
         Accumulator := Accumulator / Base;

         if Result (Position) /= 0 then
            Used := Position;
         end if;
      end loop;

      if Accumulator = 0 then
         return Result (Result'First .. Used);
      else
         raise Constraint_Error;
      end if;
   end "+";

   function "*" (Left, Right : in Number) return Number is
      Accumulator : Accumulated_Value;
      Result      : Number (Left'First + Right'First ..
                            Left'Last  + Right'Last + 1) := (others => 0);
      Used        : Integer := Integer'First;
   begin
      for L in Left'Range loop
         for R in Right'Range loop
            Accumulator := Left (L) * Right (R);

            for Position in L + R .. Result'Last loop
               exit when Accumulator = 0;

               Accumulator := Accumulator + Result (Position);
               Result (Position) := Accumulator mod Base;
               Accumulator := Accumulator / Base;
               Used := Position;
            end loop;
         end loop;
      end loop;

      return Result (Result'First .. Used);
   end "*";

   procedure Divide (Dividend  : in     Number;
                     Divisor   : in     Digit;
                     Result    :    out Number;
                     Remainder :    out Digit) is
      Accumulator : Accumulated_Value := 0;
   begin
      Result := (others => 0);

      for Position in reverse Dividend'Range loop
         Accumulator := Accumulator * Base + Dividend (Position);
         Result (Position) := Accumulator / Divisor;
         Accumulator := Accumulator mod Divisor;
      end loop;

      Remainder := Accumulator;
   end Divide;

   function Trim (Item : in Number) return Number is
   begin
      for Position in reverse Item'Range loop
         if Item (Position) /= 0 then
            return Item (Item'First .. Position);
         end if;
      end loop;

      return Zero;
   end Trim;
end Long_Multiplication;

And finally we have the requested test application:

with Ada.Text_IO;
with Long_Multiplication;

procedure Test_Long_Multiplication is
   use Ada.Text_IO, Long_Multiplication;

   N : Number := Value ("18446744073709551616");
   M : Number := N * N;
begin
   Put_Line (Image (N) & " * " & Image (N) & " = " & Image (M));
end Test_Long_Multiplication;
Output:
18446744073709551616 * 18446744073709551616 = 340282366920938463463374607431768211456

Using modular types

The following implementation uses representation of a long number by an array of 32-bit elements:

type Long_Number is array (Natural range <>) of Unsigned_32;

function "*" (Left, Right : Long_Number) return Long_Number is
   Result : Long_Number (0..Left'Length + Right'Length - 1) := (others => 0);
   Accum  : Unsigned_64;
begin
   for I in Left'Range loop
      for J in Right'Range loop
         Accum := Unsigned_64 (Left (I)) * Unsigned_64 (Right (J));
         for K in I + J..Result'Last loop
            exit when Accum = 0;
            Accum := Accum + Unsigned_64 (Result (K));
            Result (K) := Unsigned_32 (Accum and 16#FFFF_FFFF#);
            Accum := Accum / 2**32;
         end loop;
      end loop;
   end loop;
   for Index in reverse Result'Range loop -- Normalization
      if Result (Index) /= 0 then
         return Result (0..Index);
      end if;
   end loop;
   return (0 => 0);
end "*";

The task requires conversion into decimal base. For this we also need division to short number with a remainder. Here it is:

procedure Div
          (  Dividend  : in out Long_Number;
             Last      : in out Natural;
             Remainder : out Unsigned_32;
             Divisor   : Unsigned_32
          )  is
   Div   : constant Unsigned_64 := Unsigned_64 (Divisor);
   Accum : Unsigned_64 := 0;
   Size  : Natural     := 0;
begin
   for Index in reverse Dividend'First..Last loop
      Accum := Accum * 2**32 + Unsigned_64 (Dividend (Index));
      Dividend (Index) := Unsigned_32 (Accum / Div);
      if Size = 0 and then Dividend (Index) /= 0 then
         Size := Index;
      end if;
      Accum := Accum mod Div;
   end loop;
   Remainder := Unsigned_32 (Accum);
   Last := Size;
end Div;

With the above the test program:

with Ada.Strings.Unbounded;  use Ada.Strings.Unbounded;
with Ada.Text_IO;            use Ada.Text_IO;
with Interfaces;             use Interfaces;

procedure Long_Multiplication is
   -- Insert definitions above here
   procedure Put (Value : Long_Number) is
      X      : Long_Number := Value;
      Last   : Natural     := X'Last;
      Digit  : Unsigned_32;
      Result : Unbounded_String;
   begin
      loop
         Div (X, Last, Digit, 10);
         Append (Result, Character'Val (Digit + Character'Pos ('0')));
         exit when Last = 0 and then X (0) = 0;
      end loop;
      for Index in reverse 1..Length (Result) loop
         Put (Element (Result, Index));
      end loop;
   end Put;
   
   X : Long_Number := (0 => 0, 1 => 0, 2 => 1) * (0 => 0, 1 => 0, 2 => 1);
begin
   Put (X);
end Long_Multiplication;

Sample output:

340282366920938463463374607431768211456

Aime

data b, c, v;
integer d, e, i, j, s;

b = 1.argv;
b.dump(',');
v = 2.argv;
v.dump(',');

c.run(~b + ~v + 1, 0);

for (i, d in b) {
    b[i] = d - '0';
}

for (j, d of v) {
    d = v[j] - '0';

    s = 0;
    for (i, e of b) {
        s += e * d + c[i + j];
        c[i + j] = s % 10;
        s /= 10;
    }
    while (s) {
        s += c[i + j];
        c[i + j] = s % 10;
        s /= 10;
        i -= 1;
    }
}

c.delete(-1);
c.bf_drop0("");

for (i, d in c) {
    c[i] = d + '0';
}

o_form("~\n", c);

ALGOL 60

Works with: GNU MARST version 2.7
Translation of: ATS
procedure multiplyBCD (m, n, u, v, w);
  comment
    Multiply array u of length m by array v of length n,
    putting the result in array w of length (m + n). the
    numbers are stored as binary coded decimal, most
    significant digit first.
  ;
  value m, n;
  integer m, n;
  integer array u, v, w;
begin
  integer i, j, carry, t;

  for j := 0 step 1 until n - 1 do
    begin
      if v[n - 1 - j] = 0 then
        begin
          comment (optional branch);
          w[n - 1 - j] := 0
        end
      else
        begin
          carry := 0;
          for i := 0 step 1 until m - 1 do
            begin
              t := (u[m - 1 - i] * v[n - 1 - j])
                      + w[m + n - 1 - i - j] + carry;
              carry := t % 10;  comment (integer division);
              w[m + n - 1 - i - j] := t - (carry * 10)
            end;
          w[n - 1 - j] := carry
        end
    end
end;

procedure printBCD (m, u);
  value m;
  integer m;
  integer array u;
begin
  integer i, j;

  comment  Skip leading zeros;
  i := 0;
  for j := i while j < m - 1 & u[j] = 0 do
    i := i + 1;

  comment  Print the digits, and separators;
  for j := i step 1 until m - 1 do
    begin
      if j != i & ((m - j) % 3) * 3 = m - j then
        begin
          comment  Print UTF-8 for a narrow no-break space (U+202F);
          outstring (1, "\xE2\x80\xAF")
        end;
      outchar (1, "0123456789", u[j] + 1)
    end
end;

begin
  integer array u[0 : 19];
  integer array v[0 : 19];
  integer array w[0 : 39];

  u[0] := 1; u[1] := 8; u[2] := 4; u[3] := 4;
  u[4] := 6; u[5] := 7; u[6] := 4; u[7] := 4;
  u[8] := 0; u[9] := 7; u[10] := 3; u[11] := 7;
  u[12] := 0; u[13] := 9; u[14] := 5; u[15] := 5;
  u[16] := 1; u[17] := 6; u[18] := 1; u[19] := 6;

  v[0] := 1; v[1] := 8; v[2] := 4; v[3] := 4;
  v[4] := 6; v[5] := 7; v[6] := 4; v[7] := 4;
  v[8] := 0; v[9] := 7; v[10] := 3; v[11] := 7;
  v[12] := 0; v[13] := 9; v[14] := 5; v[15] := 5;
  v[16] := 1; v[17] := 6; v[18] := 1; v[19] := 6;

  multiplyBCD (20, 20, u, v, w);
  outstring (1, "u = "); printBCD (20, u); outstring (1, "\n");
  outstring (1, "v = "); printBCD (20, v); outstring (1, "\n");
  outstring (1, "u × v = "); printBCD (40, w); outstring (1, "\n")
end
Output:
$ marst long_mult_task.algol60 > algol60-code.c && cc algol60-code.c -lalgol && ./a.out
u = 18 446 744 073 709 551 616
v = 18 446 744 073 709 551 616
u × v = 340 282 366 920 938 463 463 374 607 431 768 211 456

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.

PRAGMAT precision=200 PRAGMAT
MODE INTEGER = LONG LONG INT;

LONG INT default integer width := 69;
INT width = 69+2;

INT fix w = 1, fix h = 1; # round up #

LONG LONG INT golden ratio w := ENTIER ((long long sqrt(5)-1) / 2 * LENG LENG 10 ** default integer width + fix w),
              golden ratio h := ENTIER ((long long sqrt(5)+1) / 2 * LENG LENG 10 ** default integer width + fix h);

test: (
  print((
    "The approximate golden ratios, width: ",  whole(golden ratio w,width), new line,
    "                              length: ", whole(golden ratio h,width), new line,
    "                product is exactly: ", whole(golden ratio w*golden ratio h,width*2), new line));

  INTEGER two to the power of 64 = LONG 2 ** 64;
  INTEGER neg two to the power of 64 = -(LONG 2 ** 64);
  print(("2 ** 64 * -(2 ** 64) = ", whole(two to the power of 64*neg two to the power of 64,width), new line))
)

Output:

The approximate golden ratios, width:  +618033988749894848204586834365638117720309179805762862135448622705261
                              length: +1618033988749894848204586834365638117720309179805762862135448622705261
                product is exactly:   +1000000000000000000000000000000000000000000000000000000000000000000001201173450350400438606015942314498798603569682901026716145698077078121
2 ** 64 * -(2 ** 64) =                                -340282366920938463463374607431768211456

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

# "digits" are stored in digit base ten, but 10000 & 2**n (inc hex) can be used #
INT digit base = 1000;

# if possible, then print the digit with one character #
STRING hex digit repr = "0123456789abcdefghijklmnopqrstuvwxyz"[AT 0];
INT digit base digit width = ( digit base <=  UPB hex digit repr + 1 | 1 | 1 + ENTIER log(digit base-1) );

INT next digit = -1; # reverse order so digits appear in "normal" order when printed #

PROC raise value error = ([]STRING args)VOID: 
  ( print(("Value Error: ", args, new line)); stop );

PROC raise not implemented error = ([]STRING args)VOID: 
  ( print(("Not implemented Error: ", args, new line)); stop );

PROC raise integer not implemented error = (STRING message)INTEGER: 
  ( raise not implemented error(("INTEGER ", message)); SKIP );

INT half max int = max int OVER 2;
IF digit base > half max int THEN raise value error("INTEGER addition may fail") FI;

INT sqrt max int = ENTIER sqrt(max int);
IF digit base > sqrt max int THEN raise value error("INTEGER multiplication may fail") FI;

# initialise/cast a INTEGER from a LONG LONG INT #
OP INTEGERINIT = (LONG LONG INT number)INTEGER:(
  [1 + ENTIER (SHORTEN SHORTEN long long log(ABS number) / log(digit base))]DIGIT out;
  LONG LONG INT carry := number;
  FOR digit out FROM UPB out BY next digit TO LWB out DO
    LONG LONG INT prev carry := carry;
    carry %:= digit base; # avoid MOD as it doesn't under handle -ve numbers #
    out[digit out] := SHORTEN SHORTEN (prev carry - carry * digit base)
  OD;
  out
);

# initialise/cast a INTEGER from an LONG INT #
OP INTEGERINIT = (LONG INT number)INTEGER: INTEGERINIT LENG number;

# initialise/cast a INTEGER from an INT #
OP INTEGERINIT = (INT number)INTEGER: INTEGERINIT LENG LENG number;

# remove leading zero "digits" #
OP NORMALISE = ([]DIGIT number)INTEGER: (
  INT leading zeros := LWB number - 1;
  FOR digit number FROM LWB number TO UPB number 
    WHILE number[digit number] = 0 DO leading zeros := digit number OD;
  IF leading zeros = UPB number THEN 0 ELSE number[leading zeros+1:] FI
);

#####################################################################
  Define a standard representation for the INTEGER mode.  Note: this is
  rather crude because for a large "digit base" the number is represented as
  blocks of decimals. It works nicely for powers of ten (10,100,1000,...),
  but for most larger bases (greater then 35) the repr will be a surprise.
#####################################################################
OP REPR = (DIGIT d)STRING:
    IF digit base > UPB hex digit repr THEN
      STRING out := whole(ABS d, -digit base digit width);
# Replace spaces with zeros #
      FOR digit out FROM LWB out TO UPB out DO
        IF out[digit out] = " " THEN out[digit out] := "0" FI
      OD;
      out
    ELSE # small enough to represent as ASCII (hex) characters #
      hex digit repr[ABS d]
    FI;

OP REPR = (INTEGER number)STRING:(
  STRING sep = ( digit base digit width > 1 | "," | "" );
  INT width := digit base digit width + UPB sep;
  [width * UPB number - UPB sep]CHAR out;
  INT leading zeros := LWB out - 1; 
  FOR digit TO UPB number DO
    INT start := digit * width - width + 1;
    out[start:start+digit base digit width-1] := REPR number[digit];
    IF digit base digit width /= 1 & digit /= UPB number THEN
      out[start+digit base digit width] := ","
    FI
  OD;

# eliminate leading zeros #
  FOR digit out FROM LWB out TO UPB out 
    WHILE out[digit out] = "0" OR out[digit out] = sep 
  DO leading zeros := digit out OD;

  CHAR sign = ( number[1]<0 | "-" | "+" );
# finally return the semi-normalised result #
  IF leading zeros = UPB out THEN "0" ELSE sign + out[leading zeros+1:] FI
);
################################################################
# Finally Define the required INTEGER multiplication OPerator. #
################################################################
OP * = (INTEGER a, b)INTEGER:(

# initialise out to all zeros #
  [UPB a + UPB b]INT ab; FOR place ab TO UPB ab DO ab[place ab]:=0 OD; 

  FOR place a FROM UPB a BY next digit TO LWB a DO
    DIGIT carry := 0;

# calculate each digit (whilst removing the carry) #
    FOR place b FROM UPB b BY next digit TO LWB b DO
      # n.b. result may be 2 digits #
      INT result := ab[place a + place b] + a[place a]*b[place b] + carry;
      carry := result % digit base; # avoid MOD as it doesn't under handle -ve numbers #
      ab[place a + place b] := result  - carry * digit base
    OD;
    ab[place a + LWB b + next digit] +:= carry

  OD;
  NORMALISE ab
);
# The following standard operators could (potentially) also be defined #
OP -   = (INTEGER a)INTEGER: raise integer not implemented error("monadic minus"),
  ABS  = (INTEGER a)INTEGER: raise integer not implemented error("ABS"),
  ODD  = (INTEGER a)INTEGER: raise integer not implemented error("ODD"),
  BIN  = (INTEGER a)INTEGER: raise integer not implemented error("BIN");

OP +  = (INTEGER a, b)INTEGER: raise integer not implemented error("addition"),
   -  = (INTEGER a, b)INTEGER: raise integer not implemented error("subtraction"),
   /  = (INTEGER a, b)REAL: ( VOID(raise integer not implemented error("floating point division")); SKIP),
   %  = (INTEGER a, b)INTEGER: raise integer not implemented error("fixed point division"),
   %* = (INTEGER a, b)INTEGER: raise integer not implemented error("modulo division"),
   ** = (INTEGER a, b)INTEGER: raise integer not implemented error("to the power of");

LONG INT default integer width := long long int width - 2;

INT fix w = -1177584, fix h = -3915074; # floating point error, probably GMP/hardware specific #

INTEGER golden ratio w := INTEGERINIT ENTIER ((long long sqrt(5)-1) / 2 * LENG LENG 10 ** default integer width + fix w),
        golden ratio h := INTEGERINIT ENTIER ((long long sqrt(5)+1) / 2 * LENG LENG 10 ** default integer width + fix h);

test: (
  print((
    "The approximate golden ratios, width: ",  REPR golden ratio w, new line,
    "                            length: ", REPR golden ratio h, new line,
    "                product is exactly: ", REPR (golden ratio w * golden ratio h), new line));

  INTEGER two to the power of 64 = INTEGERINIT(LONG 2 ** 64);
  INTEGER neg two to the power of 64 = INTEGERINIT(-(LONG 2 ** 64));
  print(("2 ** 64 * -(2 ** 64) = ", REPR (two to the power of 64 * neg two to the power of 64), new line))
)

Output:

The approximate golden ratios, width: +618,033,988,749,894,848,204,586,834,365,638,117,720,309,179,805,762,862,135,448,622,705,261
                            length: +1,618,033,988,749,894,848,204,586,834,365,638,117,720,309,179,805,762,862,135,448,622,705,261
                product is exactly: +1,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,001,201,173,450,350,400,438,606,015,942,314,498,798,603,569,682,901,026,716,145,698,077,078,121
2 ** 64 * -(2 ** 64) = -340,282,366,920,938,463,463,374,607,431,768,211,456

Other libraries or implementation specific extensions

As of February 2009 no open source libraries to do this task have been located.

ALGOL W

begin
    % long multiplication of large integers                                 %
    % large integers are represented by arrays of integers whose absolute   %
    % values are in 0 .. ELEMENT_MAX - 1                                    %
    % negative large integers should have negative values in all non-zero   %
    % elements                                                              %
    % the least significant digits of the large integer are in element 1    %
    integer ELEMENT_DIGITS; % number of digits in an element of a large     %
                            % integer                                       %
    integer ELEMENT_MAX;    % max absolute value of an element of a large   %
                            % integer - must be 10^( ELEMENT_DIGITS + 1 )   %
    integer ELEMENT_COUNT;  % number of elements in each large integer      %
    % implements long multiplication, c is set to a * b                     %
    %     c can be the same array as a or b                                 %
    %     n is the number of elements in the large integers a, b and c      %
    procedure longMultiply( integer array a, b, c ( * )
                          ; integer value n
                          ) ;
    begin
        % multiplies the large integer in b by the integer a, the result    %
        %     is added to c, starting from offset                           %
        %     overflow is ignored                                           %
        procedure multiplyElement( integer value a
                                 ; integer array b, c ( * )
                                 ; integer value offset, n
                                 ) ;
        begin
            integer carry, cPos;
            carry := 0;
            cPos  := offset;
            for bPos := 1 until highestNonZeroElementPosition( b, ( n + 1 ) - offset ) do begin
                integer cElement;
                cElement := c( cPos ) + ( a * b( bPos ) ) + carry;
                if abs cElement < ELEMENT_MAX then carry := 0
                else begin
                    % have digits to carry                                  %
                    carry    := cElement div ELEMENT_MAX;
                    cElement := ( abs cElement ) rem ELEMENT_MAX;
                    if carry < 0 then cElement := - cElement
                end if_no_carry_ ;
                c( cPos ) := cElement;
                cPos := cPos + 1
            end for_aPos ;
            if cPos <= n then c( cPos ) := carry
        end multiplyElement ;
        integer array mResult ( 1 :: n );
        % the result will be computed in mResult, allowing a or b to be c   %
        for rPos := 1 until n do mResult( rPos ) := 0;
        % multiply and add each element to the result                       %
        for aPos := 1 until highestNonZeroElementPosition( a, n ) do begin
            if a( aPos ) not = 0 then multiplyElement( a( aPos ), b, mResult, aPos, n )
        end for_aPos ;
        % return the result in c                                            %
        for rPos := 1 until n do c( rPos ) := mResult( rPos )
    end longMultiply ;
    % writes the decimal value of a large integer a with n elements         %
    procedure writeonLargeInteger( integer array a ( * )
                                 ; integer value n
                                 ) ;
    begin
        integer aMax;
        aMax := highestNonZeroElementPosition( a, n );
        if aMax < 1 then writeon( "0" )
        else begin
            % the large integer is non-zero                                 %
            writeon( i_w := 1, s_w := 0, a( aMax ) ); % highest element     %
            % handle the remaining elements - show leading zeros            %
            for aPos := aMax - 1 step -1 until 1 do begin
                integer v;
                integer array digits ( 1 :: ELEMENT_DIGITS );
                v := abs a( aPos );
                for dPos := ELEMENT_DIGITS step -1 until 1 do begin
                    digits( dPos ) := v rem 10;
                    v              := v div 10
                end for_dPos;
                for dPos := 1 until ELEMENT_DIGITS do writeon( i_w := 1, s_w := 0, digits( dPos ) )
            end for_aPos
        end if_aMax_lt_1_
    end writeonLargeInteger ;
    % returns the position of the highest non-zero element of the large     %
    %     integer a with n elements                                         %
    integer procedure highestNonZeroElementPosition( integer array a ( * )
                                                   ; integer value n
                                                   ) ;
    begin
        integer aMax;
        aMax := n;
        while aMax > 0 and a( aMax ) = 0 do aMax := aMax - 1;
        aMax
    end highestNonZeroElementPosition ;
    % allow each element to contain 4 decimal digits, so element by element %
    % multiplication won't overflow 32-bits                                 %
    ELEMENT_DIGITS :=     4;
    ELEMENT_MAX    := 10000;
    ELEMENT_COUNT  :=    12; % allows up to 48 digits - enough for the task %
    begin
        integer array twoTo64, twoTo128 ( 1 :: ELEMENT_COUNT );
        integer pwr;
        % construct 2^64 in twoTo64                                         %
        for tPos := 2 until ELEMENT_COUNT do twoTo64( tPos ) := 0;
        twoTo64( 1 ) := 2;
        pwr          := 1;
        while pwr < 64 do begin
            longMultiply( twoTo64, twoTo64, twoTo64, ELEMENT_COUNT );
            pwr := pwr * 2
        end while_pwr_lt_64 ;
        % construct 2^128                                                   %
        longMultiply( twoTo64, twoTo64, twoTo128, ELEMENT_COUNT );
        write( "2^128: " );
        writeonLargeInteger( twoTo128, ELEMENT_COUNT )
    end
end.
Output:
2^128: 340282366920938463463374607431768211456

APL

Works with: Dyalog APL

This function takes two digit vectors of arbitrary size.

longmul{⎕IO0
    sz/¨x y↓⌽↑⌽¨⍺⍵
    ds+(⌽⍳sz)¨(¯2×sz)[1]x∘.×y
    mlt{(1⌽⌊÷10)+10|}≡⊢ds
    0=≢mlt(\0mlt)/mlt:,0
    mlt
}
Output:
      inputlongmul63⊢,2  ⍝ construct 2*64
1 8 4 4 6 7 4 4 0 7 3 7 0 9 5 5 1 6 1 6
      longmulinput         ⍝ calculate 2*128
3 4 0 2 8 2 3 6 6 9 2 0 9 3 8 4 6 3 4 6 3 3 7 4 6 0 7 4 3 1 7 6 8 2 1 1 4 5 6

ARM Assembly

Works with: as version Raspberry Pi
or android 32 bits with application Termux
/* ARM assembly Raspberry PI  */
/*  program longmulti.s   */

 /* REMARK 1 : this program use routines in a include file 
   see task Include a file language arm assembly 
   for the routine affichageMess conversion10 
   see at end of this program the instruction include */
   
/* REMARK 2 : this program use factors unsigned to 2 power 95
              and the result is less than 2 power 159 */
   
   
/* for constantes see task include a file in arm assembly */
/************************************/
/* Constantes                       */
/************************************/
.include "../constantes.inc"
.equ BUFFERSIZE,   64

/***********************************************/
/* structures                                  */
/**********************************************/
/* Définition multi128 */
    .struct  0
multi128_N1:                      //   31-0
    .struct  multi128_N1 + 4  
multi128_N2:                      //  63-32
    .struct  multi128_N2 + 4 
multi128_N3:                      //  95-64
    .struct  multi128_N3 + 4 
multi128_N4:                      //  127-96
    .struct  multi128_N4 + 4 
multi128_N5:                      //  159-128
    .struct  multi128_N5 + 4 
multi128_end:
/*********************************/
/* Initialized data              */
/*********************************/
.data
szMessFactor:         .asciz "Factor = "
szMessResult:         .asciz "Result = "
szMessStart:          .asciz "Program 32 bits start.\n"
szCarriageReturn:     .asciz "\n"

i128test1:            .int  0,0,1,0,0    // 2 power 64
/*********************************/
/* UnInitialized data            */
/*********************************/
.bss
sZoneConv:             .skip BUFFERSIZE   // conversion buffer
i128Result1:           .skip multi128_end 

/*********************************/
/*  code section                 */
/*********************************/
.text
.global main 
main:                            @ entry of program 
    ldr r0,iAdrszMessStart
    bl affichageMess
    ldr r0,iAdri128test1         @ origin number
    ldr r1,iAdrsZoneConv
    mov r2,#BUFFERSIZE
    bl convertMultiForString     @ convert multi number to string
    mov r2,r0                    @ insert conversion in message
    mov r0,#3                    @ string number to display
    ldr r1,iAdrszMessFactor
    ldr r3,iAdrszCarriageReturn
    bl displayStrings            @ display message
                                 @ multiplication
    ldr r0,iAdri128test1         @ factor 1
    ldr r1,iAdri128test1         @ factor 2
    ldr r2,iAdri128Result1       @ result 
    bl multiplierMulti128
    ldr r0,iAdri128Result1
    ldr r1,iAdrsZoneConv
    mov r2,#BUFFERSIZE
    bl convertMultiForString     @ conversion multi to string
    mov r2,r0                    @ insert conversion in message
    mov r0,#3                    @ number string to display
    ldr r1,iAdrszMessResult
    ldr r3,iAdrszCarriageReturn
    bl displayStrings            @ display message
 
100:                             @ standard end of the program 
    mov r0, #0                   @ return code
    mov r7, #EXIT                @ request to exit program
    svc #0                       @ perform the system call
    
iAdrszCarriageReturn:        .int szCarriageReturn
iAdrsZoneConv:               .int sZoneConv
iAdri128test1:               .int i128test1
iAdri128Result1:             .int i128Result1
iAdrszMessResult:            .int szMessResult
iAdrszMessFactor:            .int szMessFactor
iAdrszMessStart:             .int szMessStart
/***************************************************/
/*  multiplication multi128 by multi128             */
/***************************************************/
// r0 contains address multi128 1
// r1 contains address multi128 2
// r2 contains address result multi128
// r0 return address result (= r2)
multiplierMulti128:
    push {r1-r10,lr}             @ save registers 
    mov r9,r0                    @ factor 1
    mov r10,r1                   @ factor 2
    mov r7,r2                    @ address result
    mov r6,#4                    @ multi128 size 
    mov r5,#0
1:
    str r5,[r7,r6,lsl #2]        @ init result
    subs r6,r6,#1
    bge 1b
    mov r5,#0                    @ indice loop 1
2:                               @ loop items factor 1
    ldr r0,[r9,r5,lsl #2]        @ load a item
    mov r4,#0
    mov r8,#0
3:                               @ loop item factor 2
    add r6,r4,r5                 @ compute result indice
 
    ldr r1,[r10,r4,lsl #2]       @ oad a item factor 2
    umull r2,r3,r1,r0            @ multiply long 32 bits 
    ldr r1,[r7,r6,lsl #2]        @ load previous item of result
    adds r1,r1,r2                @ add low part result multiplication
    movcc r2,#0                  @ high retain
    movcs r2,#1
    adds  r1,r1,r8               @ add high part precedente
    adc  r8,r3,r2                @ new high part with retenue
    str r1,[r7,r6,lsl #2]        @ store the sum in result
   
    add r4,r4,#1
    cmp r4,#3
    blt 3b                       @ and loop  2
    cmp r8,#0                    @ high part ?
    beq 4f
    add r6,r6,#1
    cmp r6,#4                    @ on last item ?
    strle r8,[r7,r6,lsl #2]      @ no store high part in next item
    ble 4f  
    adr r0,szMessErrOverflow     @ yes -> overflow
    bl affichageMess
    mov r0,#0                    @ return 0
    b 100f
4:
    add r5,r5,#1 
    cmp r5,#3
    blt 2b                       @ and loop 1
    mov r0,r7

100:
    pop {r1-r10,pc}              @ restaur registers
szMessErrOverflow:        .asciz "\033[31mOverflow !!\033[0m \n"
.align 4
/***************************************************/
/*   display multi strings                    */
/***************************************************/
/* r0  contains number strings address */
/* r1 address string1 */
/* r2 address string2 */
/* r3 address string3 */
/* other address on the stack */
/* thinck to add  number other address * 4 to add to the stack */
displayStrings:            @ INFO:  displayStrings
    push {r1-r4,fp,lr}     @ save des registres
    add fp,sp,#24          @ save paraméters address (6 registers saved * 4 bytes)
    mov r4,r0              @ save strings number
    cmp r4,#0              @ 0 string -> end
    ble 100f
    mov r0,r1              @ string 1
    bl affichageMess
    cmp r4,#1              @ number > 1
    ble 100f
    mov r0,r2
    bl affichageMess
    cmp r4,#2
    ble 100f
    mov r0,r3
    bl affichageMess
    cmp r4,#3
    ble 100f
    mov r3,#3
    sub r2,r4,#4
1:                         @ loop extract address string on stack
    ldr r0,[fp,r2,lsl #2]
    bl affichageMess
    subs r2,#1
    bge 1b
100:
    pop {r1-r4,fp,pc}
/***************************************************/
/*   conversion multi128 unsigned to string    */
/***************************************************/
// r0 contains address multi128
// r1 contains address buffer
// r2 contains buffer length
convertMultiForString:
    push {r1-r5,fp,lr}        @ save des registres
    sub sp,sp,#multi128_end   @ reserve place to stack
    mov fp,sp                 @ init address to quotient
    mov r5,r1                 @ save address buffer
    mov r3,#0                 @ init indice
1:
    ldr r4,[r0,r3,lsl #2]     @ load one part of number
    str r4,[fp,r3,lsl #2]     @ copy part on stack
    add r3,#1
    cmp r3,#5
    blt 1b
 
2:
    mov r0,#0
    strb r0,[r5,r2]            @ store final 0 in buffer
    sub r4,r2,#1               @ end number storage
3:
    mov r0,fp
    mov r1,#10
    bl calculerModuloMultiEntier @ compute modulo 10
    add r0,r0,#0x30             @ convert result to character
    strb r0,[r5,r4]             @ store character on buffer
    subs r4,r4,#1               @ 
    blt 99f                     @  buffer too low
    ldr r0,[fp,#multi128_N1]    @ test if quotient = zero
    cmp r0,#0
    bne 3b
    ldr r0,[fp,#multi128_N2]
    cmp r0,#0
    bne 3b
    ldr r0,[fp,#multi128_N3]
    cmp r0,#0
    bne 3b
    ldr r0,[fp,#multi128_N4]
    cmp r0,#0
    bne 3b
    ldr r0,[fp,#multi128_N5]
    cmp r0,#0
    bne 3b

    add r0,r5,r4                 @ return begin number in buffer
    add r0,r0,#1            
    b 100f
99:                              @ display error if buffer est toop low
    adr r0,szMessErrBuffer
    bl affichageMess
    mov r0,#-1
100:
    add sp,sp,#multi128_end      @ stack alignement
    pop {r1-r5,fp,pc}            @ restaur registers
szMessErrBuffer:        .asciz "\033[31mBuffer de conversion trop petit !!\033[0m \n"
.align 4
/***************************************************/
/*    modulo  compute   unsigned                   */
/***************************************************/
// r0 contains address multi128
// r1 contains modulo (positive)
// r0 return  modulo
// ATTENTION : le multientier origine est modifié et contient le quotient
calculerModuloMultiEntier:   @ INFO: calculerModuloMultiEntier
    push {r1-r5,lr}        @ save des registres 
    cmp r1,#0
    ble 99f
    mov r4,r1              @ save modulo
    mov r3,#4
    mov r5,r0              @ multi128 address
    ldr r0,[r5,r3,lsl #2]  @ load last part of number in low part of 64 bits
    mov r1,#0              @ init higt part 64 bits
1:
    cmp r3,#0              @ end part ?
    ble 2f
    mov r2,r4              @ modulo
    bl division32R         @ divide r0,r1 by r2 in r0,r1 and remainder in r2
    str r0,[r5,r3,lsl #2]  @ store result part low  
    sub r3,r3,#1           @ other part ?
    ldr r0,[r5,r3,lsl #2]  @ load prev part
    mov r1,r2              @ store remainder un high part of 64 bits           
    b 1b
2:
    mov r2,r4              @ modulo
    bl division32R
    str r0,[r5]            @ stockage dans le 1er chunk
    mov r0,r2              @ return remainder
    b 100f
99:
    adr r0,szMessNegatif
    bl affichageMess
    mov r0,#-1
100:                       @ fin standard de la fonction  
    pop {r1-r5,pc}         @ restaur des registres
szMessNegatif:      .asciz "\033[31mLe diviseur doit être positif !\033[0m\n"
.align 4
/***************************************************/
/*   division 64 bits number in 2 registers by 32 bits number */
/***************************************************/
/* r0 contains  dividende  low part */
/* r1 contains  dividende  high part */
/* r2 contains divisor   */
/* r0 return quotient  low part  */
/* r1 return quotient  high part  */
/* r2 return remainder                    */
division32R:
    push {r3-r7,lr}   @ save registers
    mov r6,#0         @ init high high part of remainder !!
                      @ r1 = high part of number in high part of remainder
    mov r7,r0         @ low part of number in low part of remainder
    mov r3,#0         @ init high part quotient
    mov r4,#0         @ init low part quotient
    mov r5,#32
1:                    @ begin loop
    lsl r6,#1         @ left shift high high part of remainder
    lsls r1,#1        @ left shift high part of remainder 
    orrcs r6,#1       @ add left bit  high part on high high part
    lsls r7,#1        @ left shift low part of remainder 
    orrcs r1,#1       @ add left bit  low part on high part
    lsls r4,#1        @ left shift low part quotient 
    lsl r3,#1         @ left shift high part quotient
    orrcs r3,#1       @ add left bit low part on high part
                      @ sub divisor to high part remainder
    subs r1,r2
    sbcs  r6,#0       @ sub restraint (retenue in french)
    bmi 2f            @ result negative ?
    
                      @ positive or equal
    orr r4,#1         @ right bit quotient  to 1
    b 3f
2:                    @ negative
    orr r4,#0         @ right bit quotient to 0 
    adds r1,r2        @ and restaure the remainder to precedent value
    adc  r6,#0        @ and restraint
3:
    subs r5,#1        @ decrement indice 
    bgt 1b            @ and loop
    mov r0,r4         @ low part quotient 
    mov r2,r1         @ remainder
    mov r1,r3         @ high part quotient
100:                  @ 
    pop {r3-r7,pc}    @ restaur registers
/***************************************************/
/*      ROUTINES INCLUDE                           */
/***************************************************/
.include "../affichage.inc"
Output:
Program 32 bits start.
Factor = 18446744073709551616
Result = 340282366920938463463374607431768211456

Arturo

Translation of: ATS

The program is written to operate on strings containing digits and spaces.

; The following two functions assume the 7-bit encoding is ASCII.
char2BCD: function [c] [
  return (and (to :integer c) 15) % 10
]
BCD2char: function [i] [
  return (to :char (or i 48))
]

multiplyBCD: function [u v] [
  m: size u
  n: size v
  w: array.of: (m + n) `0`

  predm: m - 1
  predn: n - 1
  predszw: (size w) - 1

  ; Long multiplication. See Algorithm 4.3.1M in Volume 2 of Knuth,
  ; ‘The Art of Computer Programming’. Here b = 10. Only the less
  ; significant nibble of a character is considered. Thus zero can be
  ; represented by either `0` or ` `, and other digits by their
  ; respective ASCII characters.
  loop 0..predn 'j [
    vj: char2BCD v\[predn - j]
    if? vj = 0 [
      set w (predn - j) `0`
    ] else [
      carry: 0
      loop 0..predm 'i [
        ui: char2BCD u\[predm - i]
        wij: char2BCD w\[predszw - (i + j)]
        t: (ui * vj) + wij + carry
        [carry digit]: divmod t 10
        set w (predszw - (i + j)) (BCD2char digit)
      ]
      set w (predn - j) (BCD2char carry)
    ]
  ]

  return join w
]

twoRaised64: "18446744073709551616"
twoRaised128: multiplyBCD twoRaised64 twoRaised64

print twoRaised128
Output:
0340282366920938463463374607431768211456

ATS

I perform both binary long multiplication (designed for efficiency) and decimal arithmetic.

The example numbers, though fine for testing decimal arithmetic, are very bad for testing binary multiplication. For example, they never require a carry. So I have added (for binary multiplication) the example 79 228 162 514 264 337 593 543 950 335 squared equals 6 277 101 735 386 680 763 835 789 423 049 210 091 073 826 769 276 946 612 225. The first number is 96 one-bits.

(* This is Algorithm 4.3.1M in Volume 2 of Knuth, ‘The Art of Computer
   Programming’. *)

#include "share/atspre_staload.hats"

#define NIL list_nil ()
#define ::  list_cons

(********************** FOR BINARY ARITHMETIC ***********************)

(* We need to choose a radix for the multiplication, small enough that
   intermediate results can be represented, but big for efficiency. To
   stay within the POSIX types, I choose 2**32 as my radix. Thus
   ‘digits’ are stored in uint32 and intermediate results are stored
   in uint64.

   A number is stored as an array of uint32, with the least
   significant uint32 first. *)

extern fn
long_multiplication             (* Multiply u and v, giving w. *)
          {m, n : int}
          (m    : size_t m,
           n    : size_t n,
           u    : &array (uint32, m),
           v    : &array (uint32, n),
           w    : &array (uint32?, m + n) >> array (uint32, m + n))
    :<!refwrt> void

%{^
#include <stdint.h>
%}
extern castfn i2u32 : int -<> uint32
extern castfn u32_2i : uint32 -<> int
extern castfn i2u64 : int -<> uint64
extern castfn u32u64 : uint32 -<> uint64
extern castfn u64u32 : uint64 -<> uint32
macdef zero32 = i2u32 0
macdef zero64 = i2u64 0
macdef one32  = i2u32 1
macdef ten32  = i2u32 10
macdef mask32 = $extval (uint32, "UINT32_C (0xFFFFFFFF)")

(* The following implementation is precisely the algorithm suggested
   by Knuth, although specialized for b=2**32 and for unsigned
   integers of precisely 32 bits. *)
implement
long_multiplication {m, n} (m, n, u, v, w) =
  let
    (* Establish that the arrays have non-negative lengths. *)
    prval () = lemma_array_param u
    prval () = lemma_array_param v

    (* Knuth initializes only part of the w array. However, if we
       initialize ALL of w now, then we will not have to deal with
       complicated array views later. *)
    val () = array_initize_elt<uint32> (w, m + n, zero32)

    (* The following function includes proof of termination. *)
    fun
    jloop {j : nat | j <= n} .<n - j>.
          (u : &array (uint32, m),
           v : &array (uint32, n),
           w : &array (uint32, m + n),
           j : size_t j)
         :<!refwrt> void =
      if j = n then
        ()
      else if v[j] = zero32 then (* This branch is optional. *)
        begin
          w[j + m] := zero32;
          jloop (u, v, w, succ j)
        end
      else
        let
          fun
          iloop {i : nat | i <= m} .<m - i>.
                (u : &array (uint32, m),
                 v : &array (uint32, n),
                 w : &array (uint32, m + n),
                 i : size_t i,
                 k : uint64)    (* carry *)
              :<!refwrt> void =
            if i = m then
              w[j + m] := u64u32 k
            else
              let
                val t = (u32u64 u[i] * u32u64 v[j])
                            + u32u64 w[i + j] + k
              in
                (* The mask here is not actually needed, if uint32
                   really is treated by the C compiler as 32 bits. *)
                w[i + j] := (u64u32 t) land mask32;

                iloop (u, v, w, succ i, t >> 32)
              end
        in
          iloop (u, v, w, i2sz 0, zero64);
          jloop (u, v, w, succ j)
        end
  in
    jloop (u, v, w, i2sz 0)
  end

fn
big_integer_iseqz               (* Is a big integer equal to zero? *)
          {m : int}
          (m : size_t m,
           u : &array (uint32, m))
    :<!ref> bool =
  let
    prval () = lemma_array_param u
    fun
    loop {n : nat | n <= m} .<n>.
         (u : &array (uint32, m),
          n : size_t n)
        :<!ref> bool =
      if n = i2sz 0 then
        true
      else if u[pred n] = zero32 then
        loop (u, pred n)
      else
        false
  in
    loop (u, m)
  end

(* To print the number in decimal, we need division by 10. So here is
   ‘short division’: Exercise 4.3.1.16 in Volume 2 of Knuth.  *)
fn
short_division
          {m : int}
          (m : size_t m,
           u : &array (uint32, m),
           v : uint32,
           q : &array (uint32?, m) >> array (uint32, m),
           r : &uint32? >> uint32)
    :<!refwrt> void =
  let
    prval () = lemma_array_param u
    val () = array_initize_elt<uint32> (q, m, zero32)
    val () = r := zero32
    fun
    loop {i1 : nat | i1 <= m} .<i1>.
         (u  : &array (uint32, m),
          q  : &array (uint32, m),
          i1 : size_t i1,
          r  : &uint32)
        :<!refwrt> void =
      if i1 <> i2sz 0 then
        let
          val i = pred i1
          val tmp = (u32u64 r << 32) lor (u32u64 u[i])
          val tmp_q = tmp / u32u64 v and tmp_r = tmp mod (u32u64 v)
        in
          q[i] := u64u32 tmp_q;
          r := u64u32 tmp_r;
          loop (u, q, i, r)
        end
  in
    loop (u, q, m, r)
  end

fn
fprint_big_integer
          {m : int}
          (f : FILEref,
           m : size_t m,
           u : &array (uint32, m))
    : void =
  let
    fun
    loop1 (v    : &array (uint32, m),
           q    : &array (uint32, m),
           lst  : List0 char,
           i    : uint)
        : List0 char =
      let
        var r : uint32
        val () = short_division (m, v, ten32, q, r)
        val r = g1ofg0 (u32_2i r)
        val () = assertloc ((0 <= r) * (r <= 9))
        val digit = int2digit r
      in
        if big_integer_iseqz (m, q) then
          digit :: lst
        else if i = 2U then
          (* Insert UTF-8 for narrow no-break space U+202F *)
          loop1 (q, v, '\xE2' :: '\x80' :: '\xAF' :: digit :: lst, 0U)
        else
          loop1 (q, v, digit :: lst, succ i)
      end
    fun
    loop2 {n   : nat} .<n>.
          (lst : list (char, n))
        : void =
      case+ lst of
      | NIL => ()
      | hd :: tl => (fprint! (f, hd); loop2 tl)
  in
    if big_integer_iseqz (m, u) then
      fprint! (f, "0")
    else
      let
        val @(pf, pfgc | p) = array_ptr_alloc<uint32> m
        val @(qf, qfgc | q) = array_ptr_alloc<uint32> m
        val () = array_copy<uint32> (!p, u, m)
        val () = array_initize_elt<uint32> (!q, m, zero32)
        val () = loop2 (loop1 (!p, !q, NIL, 0U))
        val () = array_ptr_free (pf, pfgc | p)
        val () = array_ptr_free (qf, qfgc | q)
      in
      end
  end

fn
example_binary (f : FILEref) : void =
  let
    var u = @[uint32][3] (zero32, zero32, one32)
    var v = @[uint32][3] (zero32, zero32, one32)
    var w : @[uint32][6]
  in
    long_multiplication (i2sz 3, i2sz 3, u, v, w);
    fprint! (f, "\nBinary long multiplication (b = 2³²)\n\n");
    fprint! (f, "u = ");
    fprint_big_integer (f, i2sz 3, u);
    fprint! (f, "\nv = ");
    fprint_big_integer (f, i2sz 3, v);
    fprint! (f, "\nu × v = ");
    fprint_big_integer (f, i2sz 6, w);
    fprint! (f, "\n")
  end

fn
test_binary (f : FILEref) : void =
  let
    var u = @[uint32][3] (mask32, mask32, mask32)
    var v = @[uint32][3] (mask32, mask32, mask32)
    var w : @[uint32][6]
  in
    long_multiplication (i2sz 3, i2sz 3, u, v, w);
    fprint! (f, "\nThe example numbers specified in the task\n",
                "are actually VERY bad for testing binary\n",
                "multiplication, because they never need a carry.\n",
                "So here is a multiplication full of carries,\n",
                "with b = 2³²\n\n");
    fprint! (f, "u = ");
    fprint_big_integer (f, i2sz 3, u);
    fprint! (f, "\nv = ");
    fprint_big_integer (f, i2sz 3, v);
    fprint! (f, "\nu × v = ");
    fprint_big_integer (f, i2sz 6, w);
    fprint! (f, "\n")
  end

(************** FOR BINARY CODED DECIMAL ARITHMETIC *****************)

(* The following will operate on arrays of BCD digits, with the most
   significant digit first. Only the least four bits of a byte will be
   considered. This has at least two benefits: any ASCII digit is
   treated as its BCD equivalent, and SPACE is treated as zero. *)

extern fn
bcd_multiplication              (* Multiply u and v, giving w. *)
          {m, n : int}
          (m    : size_t m,
           n    : size_t n,
           u    : &array (char, m),
           v    : &array (char, n),
           w    : &array (char?, m + n) >> array (char, m + n))
    :<!refwrt> void

fn {}
char2bcd (c : char) :<> intBtwe (0, 9) =
  let
    val c = char2uchar1 (g1ofg0 c)
    val i = g1uint_of_uchar1 c
    val i = i mod 16U
    val i = i mod 10U           (* Guarantees the digit be BCD. *)
  in
    u2i i
  end

extern castfn bcd2char (i : intBtwe (0, 9)) :<> char

(* The following implementation is precisely the algorithm suggested
   by Knuth, specialized for b=10. *)
implement
bcd_multiplication {m, n} (m, n, u, v, w) =
  let
    (* Establish that the arrays have non-negative lengths. *)
    prval () = lemma_array_param u
    prval () = lemma_array_param v

    (* Knuth initializes only part of the w array. However, if we
       initialize ALL of w now, then we will not have to deal with
       complicated array views later. *)
    val () = array_initize_elt<char> (w, m + n, '\0')

    (* The following function includes proof of termination. *)
    fun
    jloop {j : nat | j <= n} .<n - j>.
          (u : &array (char, m),
           v : &array (char, n),
           w : &array (char, m + n),
           j : size_t j)
         :<!refwrt> void =
      if j = n then
        ()
      else if char2bcd v[pred n - j] = 0 then (* Optional branch. *)
        begin
          w[pred n - j] := '\0';
          jloop (u, v, w, succ j)
        end
      else
        let
          fun
          iloop {i : nat | i <= m} .<m - i>.
                (u : &array (char, m),
                 v : &array (char, n),
                 w : &array (char, m + n),
                 i : size_t i,
                 k : intBtwe (0, 9)) (* carry *)
              :<!refwrt> void =
            if i = m then
              w[pred n - j] := bcd2char k
            else
              let
                val ui = char2bcd u[pred m - i]
                and vj = char2bcd v[pred n - j]
                and wij = char2bcd w[pred (m + n) - (i + j)]

                val t = (ui * vj) + wij + k

                (* This will prove that 0 <= t *)
                prval [ui : int] EQINT () = eqint_make_gint ui
                prval [vj : int] EQINT () = eqint_make_gint vj
                prval [t : int] EQINT () = eqint_make_gint t
                prval () = mul_gte_gte_gte {ui, vj} ()
                prval () = prop_verify {0 <= t} ()

                (* But I do not feel like proving that t / 10 <= 9. *)
                val t_div_10 = t \ndiv 10 and t_mod_10 = t \nmod 10
                val () = $effmask_exn assertloc (t_div_10 <= 9)
              in
                w[pred (m + n) - (i + j)] := bcd2char t_mod_10;
                iloop (u, v, w, succ i, t_div_10)
              end
        in
          iloop (u, v, w, i2sz 0, 0);
          jloop (u, v, w, succ j)
        end
  in
    jloop (u, v, w, i2sz 0)
  end

fn
fprint_bcd {m : int}
           (f : FILEref,
            m : size_t m,
            u : &array (char, m))
    : void =
  let
    prval () = lemma_array_param u

    fun
    skip_zeros {i : nat | i <= m} .<m - i>.
               (u : &array (char, m),
                i : size_t i)
        :<!ref> [i : nat | i <= m] size_t i =
      if i = m then
        i
      else if char2bcd u[i] = 0 then
        skip_zeros (u, succ i)
      else
        i

    val [i : int] i = skip_zeros (u, i2sz 0)

    fun
    loop {j : int | i <= j; j <= m} .<m - j>.
         (u : &array (char, m),
          j : size_t j)
        : void =
      if j <> m then
        begin
          if j <> i && (m - j) mod (i2sz 3) = i2sz 0 then
            (* Print UTF-8 for narrow no-break space U+202F *)
            fprint! (f, "\xE2\x80\xAF");
          fprint! (f, int2digit (char2bcd u[j]));
          loop (u, succ j)
        end
  in
    if i = m then
      fprint! (f, "0")
    else
      loop (u, i)
  end

fn
string2bcd {n : int}
           (s : string n)
    : [p : agz]
      @(array_v (char, p, n), mfree_gc_v p | ptr p) =
  let
    val n = strlen s
    val @(pf, pfgc | p) = array_ptr_alloc<char> n
    implement
    array_initize$init<char> (i, x) =
      let
        val i = g1ofg0 i
        prval () = lemma_g1uint_param i
        val () = assertloc (i < n)
      in
        x := s[i]
      end
    val () = array_initize<char> (!p, n)
  in
    @(pf, pfgc | p)
  end

fn
example_bcd (f : FILEref) : void =
  let
    val s = g1ofg0 "18446744073709551616"

    val m = strlen s

    val @(pf_u, pfgc_u | p_u) = string2bcd s
    val @(pf_v, pfgc_v | p_v) = string2bcd s
    val @(pf_w, pfgc_w | p_w) = array_ptr_alloc<char> (m + m)
    macdef u = !p_u
    macdef v = !p_v
    macdef w = !p_w
  in
    bcd_multiplication (m, m, u, v, w);
    fprint! (f, "\nDecimal long multiplication (b = 10)\n\n");
    fprint! (f, "u = ");
    fprint_bcd (f, m, u);
    fprint! (f, "\nv = ");
    fprint_bcd (f, m, v);
    fprint! (f, "\nu × v = ");
    fprint_bcd (f, m + m, w);
    fprint! (f, "\n");
    array_ptr_free (pf_u, pfgc_u | p_u);
    array_ptr_free (pf_v, pfgc_v | p_v);
    array_ptr_free (pf_w, pfgc_w | p_w)
  end

(********************************************************************)

implement
main () =
  begin
    example_binary (stdout_ref);
    println! ();
    example_bcd (stdout_ref);
    println! ();
    test_binary (stdout_ref);
    println! ();
    0
  end
Output:
$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_LIBC long_mult_task.dats && ./a.out

Binary long multiplication (b = 2³²)

u = 18 446 744 073 709 551 616
v = 18 446 744 073 709 551 616
u × v = 340 282 366 920 938 463 463 374 607 431 768 211 456


Decimal long multiplication (b = 10)

u = 18 446 744 073 709 551 616
v = 18 446 744 073 709 551 616
u × v = 340 282 366 920 938 463 463 374 607 431 768 211 456


The example numbers specified in the task
are actually VERY bad for testing binary
multiplication, because they never need a carry.
So here is a multiplication full of carries,
with b = 2³²

u = 79 228 162 514 264 337 593 543 950 335
v = 79 228 162 514 264 337 593 543 950 335
u × v = 6 277 101 735 386 680 763 835 789 423 049 210 091 073 826 769 276 946 612 225

AutoHotkey

ahk discussion

MsgBox % x := mul(256,256)
MsgBox % x := mul(x,x)
MsgBox % x := mul(x,x) ; 18446744073709551616
MsgBox % x := mul(x,x) ; 340282366920938463463374607431768211456

mul(b,c) { ; <- b*c
   VarSetCapacity(a, n:=StrLen(b)+StrLen(c), 48), NumPut(0,a,n,"char")
   Loop % StrLen(c) {
      i := StrLen(c)+1-A_Index, cy := 0
      Loop % StrLen(b) {
         j := StrLen(b)+1-A_Index,
         t := SubStr(a,i+j,1) + SubStr(b,j,1) * SubStr(c,i,1) + cy
         cy := t // 10
         NumPut(mod(t,10)+48,a,i+j-1,"char")
      }
      NumPut(cy+48,a,i+j-2,"char")
   }
   Return cy ? a : SubStr(a,2)
}

AWK

Works with: gawk version 3.1.7
Works with: nawk version 20100523
Translation of: Tcl
BEGIN {
    DEBUG = 0
    n = 2^64
    nn = sprintf("%.0f", n)
    printf "2^64 * 2^64 = %.0f\n", multiply(nn, nn)
    printf "2^64 * 2^64 = %.0f\n", n*n
    exit
}

function multiply(x, y,     len_x,len_y,ax,ay,j,m,c,i,k,d,v,res,mul,result) {
    len_x = split_reverse(x, ax)
    len_y = split_reverse(y, ay)
    print_array(ax)
    print_array(ay)
    for (j=1; j<=len_y; j++) {
        m = ay[j]
        c = 0
        i = j - 1
        for (k=1; k<=len_x; k++) {
            d = ax[k]
            i++
            v = res[i]
            if (v == "") {
                append_array(res, 0)
                v = 0
            }
            mul = v + c + d*m
            c = int(mul / 10)
            v = mul % 10
            res[i] = v
        }
        append_array(res, c)
    }
    print_array(res)
    result = reverse_join(res)
    sub(/^0+/, "", result)
    return result
}

function split_reverse(x, a,        a_x) {
    split(x, a_x, "")
    return reverse_array(a_x, a)
}

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

function length_array(a,        len,i) {
    len = 0
    for (i in a) len++
    return len
}

function append_array(a, value,     len) {
    len = length_array(a)
    a[++len] = value
}

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

function print_array(a,         len,i) {
    if (DEBUG) {
        len = length_array(a)
        print "length=" len
        for (i=1; i<=len; i++) {
            printf("%s ", i%10)
        }
        print ""
        for (i=1; i<=len; i++) {
            #print i " " a[i]
            printf("%s ", a[i])
        }
        print ""
        print "===="
    }
}

outputs:

2^64 * 2^64 = 340282366920938463463374607431768211456
2^64 * 2^64 = 340282366920938463463374607431768211456

BASIC

Works with: QBasic

Version 1

'PROGRAM : BIG MULTIPLICATION VER #1
'LRCVS 01.01.2010
'THIS PROGRAM SIMPLY MAKES A MULTIPLICATION
'WITH ALL THE PARTIAL PRODUCTS.
'............................................................

DECLARE SUB A.INICIO (A$, B$)
DECLARE SUB B.STORE (CAD$, N$)
DECLARE SUB C.PIZARRA ()
DECLARE SUB D.ENCABEZADOS (A$, B$)
DECLARE SUB E.MULTIPLICACION (A$, B$)
DECLARE SUB G.SUMA ()
DECLARE FUNCTION F.INVCAD$ (CAD$)

RANDOMIZE TIMER
CALL A.INICIO(A$, B$)
CALL B.STORE(A$, "A")
CALL B.STORE(B$, "B")
CALL C.PIZARRA
CALL D.ENCABEZADOS(A$, B$)
CALL E.MULTIPLICACION(A$, B$)
CALL G.SUMA

SUB A.INICIO (A$, B$)
    CLS
'Note: Number of digits > 1000
    INPUT "NUMBER OF DIGITS  "; S
    CLS
    A$ = ""
    B$ = ""
    FOR N = 1 TO S
        A$ = A$ + LTRIM$(STR$(INT(RND * 9)))
    NEXT N
    FOR N = 1 TO S
        B$ = B$ + LTRIM$(STR$(INT(RND * 9)))
    NEXT N
END SUB

SUB B.STORE (CAD$, N$)
    OPEN "O", #1, N$
    FOR M = LEN(CAD$) TO 1 STEP -1
        WRITE #1, MID$(CAD$, M, 1)
    NEXT M
    CLOSE (1)
END SUB

SUB C.PIZARRA
    OPEN "A", #3, "R"
    WRITE #3, ""
    CLOSE (3)
    KILL "R"
END SUB

SUB D.ENCABEZADOS (A$, B$)
    LT = LEN(A$) + LEN(B$) + 1
    L$ = STRING$(LT, " ")
    OPEN "A", #3, "R"
    MID$(L$, LT - LEN(A$) + 1) = A$
    WRITE #3, L$
    CLOSE (3)
    L$ = STRING$(LT, " ")
    OPEN "A", #3, "R"
    MID$(L$, LT - LEN(B$) - 1) = "X " + B$
    WRITE #3, L$
    CLOSE (3)
END SUB

SUB E.MULTIPLICACION (A$, B$)
    LT = LEN(A$) + LEN(B$) + 1
    L$ = STRING$(LT, " ")
    C$ = ""
    D$ = ""
    E$ = ""
    CT1 = 1
    ACUM = 0
    OPEN "I", #2, "B"
    WHILE EOF(2) <> -1
        INPUT #2, B$
        OPEN "I", #1, "A"
        WHILE EOF(1) <> -1
            INPUT #1, A$
            RP = (VAL(A$) * VAL(B$)) + ACUM
            C$ = LTRIM$(STR$(RP))
            IF EOF(1) <> -1 THEN D$ = D$ + RIGHT$(C$, 1)
            IF EOF(1) = -1 THEN D$ = D$ + F.INVCAD$(C$)
            E$ = LEFT$(C$, LEN(C$) - 1)
            ACUM = VAL(E$)
        WEND
        CLOSE (1)
        MID$(L$, LT - CT1 - LEN(D$) + 2) = F.INVCAD$(D$)
        OPEN "A", #3, "R"
        WRITE #3, L$
        CLOSE (3)
        L$ = STRING$(LT, " ")
        ACUM = 0
        C$ = ""
        D$ = ""
        E$ = ""
        CT1 = CT1 + 1
    WEND
    CLOSE (2)
END SUB

FUNCTION F.INVCAD$ (CAD$)
    LCAD = LEN(CAD$)
    CADTEM$ = ""
    FOR CAD = LCAD TO 1 STEP -1
        CADTEM$ = CADTEM$ + MID$(CAD$, CAD, 1)
    NEXT CAD
    F.INVCAD$ = CADTEM$
END FUNCTION

SUB G.SUMA
    CF = 0
    OPEN "I", #3, "R"
    WHILE EOF(3) <> -1
        INPUT #3, R$
        CF = CF + 1
        AN = LEN(R$)
    WEND
    CF = CF - 2
    CLOSE (3)
    W$ = ""
    ST = 0
    ACUS = 0
    FOR P = 1 TO AN
        K = 0
        OPEN "I", #3, "R"
        WHILE EOF(3) <> -1
            INPUT #3, R$
            K = K + 1
            IF K > 2 THEN ST = ST + VAL(MID$(R$, AN - P + 1, 1))
            IF K > 2 THEN M$ = LTRIM$(STR$(ST + ACUS))
        WEND
        'COLOR 10: LOCATE CF + 3, AN - P + 1: PRINT RIGHT$(M$, 1); : COLOR 7
        W$ = W$ + RIGHT$(M$, 1)
        ACUS = VAL(LEFT$(M$, LEN(M$) - 1))
        CLOSE (3)
        ST = 0
    NEXT P

    OPEN "A", #3, "R"
    WRITE #3, " " + RIGHT$(F.INVCAD(W$), AN - 1)
    CLOSE (3)
    CLS
    PRINT "THE SOLUTION IN THE FILE: R"
END SUB

Version 2

'PROGRAM: BIG MULTIPLICATION VER # 2
'LRCVS 01/01/2010
'THIS PROGRAM SIMPLY MAKES A BIG MULTIPLICATION
'WITHOUT THE PARTIAL PRODUCTS.
'HERE SEE ONLY THE SOLUTION.
'...............................................................
CLS
PRINT "WAIT"

NA = 2000 'NUMBER OF ELEMENTS OF THE MULTIPLY.
NB = 2000  'NUMBER OF ELEMENTS OF THE MULTIPLIER.
'Solution = 4000 Exacts digits

'......................................................
OPEN "X" + ".MLT" FOR BINARY AS #1
CLOSE (1)
KILL "*.MLT"
'.....................................................
'CREATING THE MULTIPLY  >>> A
'CREATING THE MULTIPLIER >>> B
FOR N = 1 TO 2
IF N = 1 THEN F$ = "A" + ".MLT": NN = NA
IF N = 2 THEN F$ = "B" + ".MLT": NN = NB
    OPEN F$ FOR BINARY AS #1
        FOR N2 = 1 TO NN
            RANDOMIZE TIMER
            X$ = LTRIM$(STR$(INT(RND * 10)))
            SEEK #1, N2: PUT #1, N2, X$
        NEXT N2
    SEEK #1, N2
    CLOSE (1)
NEXT N
'.....................................................
OPEN "A" + ".MLT" FOR BINARY AS #1
FOR K = 0 TO 9
NUM$ = "": Z$ = "": ACU = 0: GG = NA
C$ = LTRIM$(STR$(K))
    OPEN C$ + ".MLT" FOR BINARY AS #2
        'OPEN "A" + ".MLT" FOR BINARY AS #1
            FOR N = 1 TO NA
                SEEK #1, GG: GET #1, GG, X$
                NUM$ = X$
                Z$ = LTRIM$(STR$(ACU + (VAL(X$) * VAL(C$))))
                L = LEN(Z$)
                ACU = 0
                IF L = 1 THEN NUM$ = Z$: PUT #2, N, NUM$
                IF L > 1 THEN ACU = VAL(LEFT$(Z$, LEN(Z$) - 1)): NUM$ = RIGHT$(Z$, 1): PUT #2, N, NUM$
                SEEK #2, N: PUT #2, N, NUM$
                GG = GG - 1
            NEXT N
        IF L > 1 THEN ACU = VAL(LEFT$(Z$, LEN(Z$) - 1)): NUM$ = LTRIM$(STR$(ACU)): XX$ = XX$ + NUM$: PUT #2, N, NUM$
        'CLOSE (1)
    CLOSE (2)
NEXT K
CLOSE (1)
'......................................................
ACU = 0
LT5 = 1
LT6 = LT5
OPEN "B" + ".MLT" FOR BINARY AS #1
    OPEN "D" + ".MLT" FOR BINARY AS #3
        FOR JB = NB TO 1 STEP -1
        SEEK #1, JB
        GET #1, JB, X$

            OPEN X$ + ".MLT" FOR BINARY AS #2: LF = LOF(2): CLOSE (2)

            OPEN X$ + ".MLT" FOR BINARY AS #2
                FOR KB = 1 TO LF
                    SEEK #2, KB
                    GET #2, , NUM$
                    SEEK #3, LT5
                    GET #3, LT5, PR$
                    T$ = ""
                    T$ = LTRIM$(STR$(ACU + VAL(NUM$) + VAL(PR$)))
                    PR$ = RIGHT$(T$, 1)
                    ACU = 0
                    IF LEN(T$) > 1 THEN ACU = VAL(LEFT$(T$, LEN(T$) - 1))
                    SEEK #3, LT5: PUT #3, LT5, PR$
                    LT5 = LT5 + 1
                NEXT KB
                IF ACU <> 0 THEN PR$ = LTRIM$(STR$(ACU)): PUT #3, LT5, PR$
            CLOSE (2)
        LT6 = LT6 + 1
        LT5 = LT6
        ACU = 0
        NEXT JB
    CLOSE (3)
CLOSE (1)
OPEN "D" + ".MLT" FOR BINARY AS #3: LD = LOF(3): CLOSE (3)
ER = 1
OPEN "D" + ".MLT" FOR BINARY AS #3
    OPEN "R" + ".MLT" FOR BINARY AS #4
        FOR N = LD TO 1 STEP -1
            SEEK #3, N: GET #3, N, PR$
            SEEK #4, ER: PUT #4, ER, PR$
            ER = ER + 1
        NEXT N
    CLOSE (4)
CLOSE (3)
KILL "D.MLT"
FOR N = 0 TO 9
    C$ = LTRIM$(STR$(N))
    KILL C$ + ".MLT"
NEXT N
PRINT "END"
PRINT "THE SOLUTION IN THE FILE: R.MLT"

Applesoft BASIC

 100 A$ = "18446744073709551616"
 110 B$ = A$
 120 GOSUB 400
 130 PRINT E$
 140 END

 400  REM MULTIPLY A$ * B$
 410 C$ = "":D$ = "0"
 420  FOR I =  LEN (B$) TO 1 STEP  - 1
 430 C = 0:B =  VAL ( MID$ (B$,I,1))
 440  FOR J =  LEN (A$) TO 1 STEP  - 1
 450 V = B *  VAL ( MID$ (A$,J,1)) + C
 460 C =  INT (V / 10):V = V - C * 10
 470 C$ =  STR$ (V) + C$
 480  NEXT J
 490  IF C THEN C$ =  STR$ (C) + C$
 510  GOSUB 600"ADD C$ + D$
 520 D$ = E$:C$ = "0":J =  LEN (B$) - I
 530  IF J THEN J = J - 1:C$ = C$ + "0": GOTO 530
 550  NEXT I
 560  RETURN

 600  REM ADD C$ + D$
 610 E =  LEN (D$):E$ = "":C = 0
 620  FOR J =  LEN (C$) TO 1 STEP  - 1
 630  IF E THEN D =  VAL ( MID$ (D$,E,1))
 640 V =  VAL ( MID$ (C$,J,1)) + D + C
 650 C = V > 9:V = V - 10 * C
 660 E$ =  STR$ (V) + E$
 670  IF E THEN E = E - 1:D = 0
 680  NEXT J
 700  IF E THEN V =  VAL ( MID$ (D$,E,1)) + C:C = V > 9:V = V - 10 * C:E$ =  STR$ (V) + E$:E = E - 1: GOTO 700
 720  RETURN

BASIC256

Translation of: Liberty BASIC
print "2^64"
a$ = "1"
for i = 1 to 64
	a$ = multByD$(a$, 2)
next
print a$
print "(check with native BASIC-256)"
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 BASIC-256)"
print 2^64*2^64
print "(looks OK)"
end

function max(a, b)
	if a > b then
		return a
	else
		return b
	end if
end function

function longMult$(a$, b$)
	signA = 1
	if left(a$,1) = "-" then
		a$ = mid(a$,2,1)
		signA = -1
	end if
	signB = 1
	if left(b$,1) = "-" then
		b$ = mid(b$,2,1)
		signB = -1
	end if

	c$ = ""
	t$ = ""
	shift$ = ""
	for i = length(a$) to 1 step -1
		d = fromradix((mid(a$,i,1)),10)
		t$ = multByD$(b$, d)
		c$ = addLong$(c$, t$+shift$)
		shift$ += "0"
	next
	if signA * signB < 0 then c$ = "-" + c$
	return c$
end function

function multByD$(a$, d)
	#multiply a$ by digit d
	c$ = ""
	carry = 0
	for i = length(a$) to 1 step -1
		a = fromradix((mid(a$,i,1)),10)
		c = a * d + carry
		carry = int(c/10)
		c = c mod 10
		c$ = string(c) + c$
	next
	if carry > 0 then c$ = string(carry) + c$
	return c$
end function

function addLong$(a$, b$)
	#add a$ + b$, for now only positive
	l = max(length(a$), length(b$))
	a$ = pad$(a$,l)
	b$ = pad$(b$,l)
	c$ = "" #result
	carry = 0
	for i = l to 1 step -1
		a = fromradix((mid(a$,i,1)),10)
		b = fromradix((mid(b$,i,1)),10)
		c = a + b + carry
		carry = int(c/10)
		c = c mod 10
		c$ = string(c) + c$
	next
	if carry > 0 then c$ = string(carry) + c$
	return c$
end function

function pad$(a$,n)  #pad$ from right with 0 to length n
	pad$ = a$
	while length(pad$) < n
		pad$ = "0" + pad$
	end while
end function
Output:
2^64
18446744073709551616
(check with native BASIC-256)
1.84467440737e+19
(looks OK)

2^64*2^64
340282366920938463463374607431768211456
(check with native BASIC-256)
3.40282366921e+38
(looks OK)

Batch File

Based on the JavaScript iterative code.

::Long Multiplication Task from Rosetta Code
::Batch File Implementation

@echo off
call :longmul 18446744073709551616 18446744073709551616 answer
echo(%answer%
exit /b 0

rem The Hellish Procedure
rem Syntax: call :longmul <n1> <n2> <variable to store product>
:longmul
    setlocal enabledelayedexpansion

    rem Define variables
    set "num1=%1"
    set "num2=%2"
    set "limit1=-1"
    set "limit2=-1"
    set "length=0"
    set "prod="

    rem Reverse the digits of each factor
    for %%A in (1,2) do (
        for /l %%B in (0,1,9) do set "num%%A=!num%%A:%%B=%%B !"
        for %%C in (!num%%A!) do ( set /a limit%%A+=1 & set "rev%%A=%%C!rev%%A!" )
    )

    rem Do the multiplication
    for /l %%A in (0,1,%limit1%) do (
        for /l %%B in (0,1,%limit2%) do (
            set /a iter=%%A+%%B
            set /a iternext=iter+1
            set /a iternext2=iter+2

            set /a prev=digit!iter!
            set /a digit!iter!=!rev1:~%%A,1!*!rev2:~%%B,1!

            rem The next line updates the length of "digits"
            if !iternext! gtr !length! set length=!iternext!
            if !iter! lss !length! set /a digit!iter!+=prev

            set /a currdigit=digit!iter!
            if !currDigit! gtr 9 (
                set /a prev=digit!iternext!
                set /a digit!iternext!=currdigit/10
                set /a digit!iter!=currdigit%%10

                rem The next line updates the length of "digits"
                if !iternext2! gtr !length! set length=!iternext2!
                if !iternext! lss !length! set /a digit!iternext!+=prev
            )
        )
    )

    rem Finalize product reversing the digits
    for /l %%F in (0,1,%length%) do set "prod=!digit%%F!!prod!"
    endlocal & set "%3=%prod%"
goto :eof
Output:
340282366920938463463374607431768211456

BBC BASIC

Library method:

      INSTALL @lib$+"BB4WMAPMLIB"
      MAPM_DllPath$ = @lib$+"BB4WMAPM.DLL"
      PROCMAPM_Init
      
      twoto64$ = "18446744073709551616"
      PRINT "2^64 * 2^64 = " ; FNMAPM_Multiply(twoto64$, twoto64$)

Explicit method:

      twoto64$ = "18446744073709551616"
      PRINT "2^64 * 2^64 = " ; FNlongmult(twoto64$, twoto64$)
      END
      
      DEF FNlongmult(num1$, num2$)
      LOCAL C%, I%, J%, S%, num1&(), num2&(), num3&()
      S% = LEN(num1$)+LEN(num2$)
      DIM num1&(S%), num2&(S%), num3&(S%)
      IF LEN(num1$) > LEN(num2$) SWAP num1$,num2$
      $$^num1&(1) = num1$
      num1&() AND= 15
      FOR I% = LEN(num1$) TO 1 STEP -1
        $$^num2&(I%) = num2$
        num2&() AND= 15
        num3&() += num2&() * num1&(I%)
        IF I% MOD 3 = 1 THEN
          C% = 0
          FOR J% = S%-1 TO I%-1 STEP -1
            C% += num3&(J%)
            num3&(J%) = C% MOD 10
            C% DIV= 10
          NEXT
        ENDIF
      NEXT I%
      num3&() += &30
      num3&(S%) = 0
      IF num3&(0) = &30 THEN = $$^num3&(1)
      = $$^num3&(0)

C

Doing it as if by hand.

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

C#

Works with: C# version 4+
If you strip out the BigInteger checking, it will work with lesser versions.

This uses the decimal type, (which has a MaxValue of 79,228,162,514,264,337,593,543,950,335). By limiting it to 10^28, it allows 28 decimal digits for the hi part, and 28 decimal digits for the lo part, 56 decimal digits total. A side computation of BigInteger assures that the results are accurate.
using System;
using static System.Console;
using BI = System.Numerics.BigInteger;

class Program {

  static decimal mx = 1E28M, hm = 1E14M, a;

  // allows for 56 digit representation, using 28 decimal digits from each decimal
  struct bi { public decimal hi, lo; }

  // sets up for squaring process
  static bi set4sq(decimal a) { bi r; r.hi = Math.Floor(a / hm); r.lo = a % hm; return r; }

  // outputs bi structure as string, optionally inserting commas
  static string toStr(bi a, bool comma = false) {
    string r = a.hi == 0 ? string.Format("{0:0}", a.lo) :
                           string.Format("{0:0}{1:" + new string('0', 28) + "}", a.hi, a.lo);
    if (!comma) return r;  string rc = "";
    for (int i = r.Length - 3; i > 0; i -= 3) rc = "," + r.Substring(i, 3) + rc;
    return r.Substring(0, ((r.Length + 2) % 3) + 1) + rc; }

  // needed because Math.Pow() returns a double
  static decimal Pow_dec(decimal bas, uint exp) {
    if (exp == 0) return 1M; decimal tmp = Pow_dec(bas, exp >> 1); tmp *= tmp;
    if ((exp & 1) == 0) return tmp; return tmp * bas; }

  static void Main(string[] args) {
    for (uint p = 64; p < 95; p += 30) {        // show prescribed output and maximum power of 2 output
      bi x = set4sq(a = Pow_dec(2M, p)), y;     // setup for squaring process
      WriteLine("The square of (2^{0}):                    {1,38:n0}", p, a); BI BS = BI.Pow((BI)a, 2);
      y.lo = x.lo * x.lo; y.hi = x.hi * x.hi;   // square lo and hi parts
      a = x.hi * x.lo * 2M;                     // calculate midterm
      y.hi += Math.Floor(a / hm);               // increment hi part w/ high part of midterm
      y.lo += (a % hm) * hm;                    // increment lo part w/ low part of midterm
      while (y.lo > mx) { y.lo -= mx; y.hi++; } // check for overflow, adjust both parts as needed
      WriteLine(" is {0,75} (which {1} match the BigInteger computation)\n", toStr(y, true),
          BS.ToString() == toStr(y) ? "does" : "fails to"); } }

}
Output:
The square of (2^64):             18,446,744,073,709,551,616
 is                         340,282,366,920,938,463,463,374,607,431,768,211,456 (which does match the BigInteger computation)

The square of (2^94): 19,807,040,628,566,084,398,385,987,584
 is 392,318,858,461,667,547,739,736,838,950,479,151,006,397,215,279,002,157,056 (which does match the BigInteger computation)

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" );
}
//--------------------------------------------------------------------------------------------------
Output:
340282366920938463463374607431768211456

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


Version 2

#include <iostream>
#include <vector>
using namespace std;

typedef unsigned long native_t;

struct ZPlus_	// unsigned int, represented as digits base 10
{
	vector<native_t> digits_;	// least significant first; value is sum(digits_[i] * 10^i)

	ZPlus_(native_t n) : digits_(1, n)
	{
		while(Sweep());
	}

	bool Sweep()	// clean up digits so they are in [0,9]
	{
		bool changed = false;
		int carry = 0;
		for (auto pd = digits_.begin(); pd != digits_.end(); ++pd)
		{
			*pd += carry;
			carry = *pd / 10;
			*pd -= 10 * carry;
			changed = changed || carry > 0;
		}
		if (carry)
			digits_.push_back(carry);
		return changed || carry > 9;
	}
};

ZPlus_ operator*(const ZPlus_& lhs, const ZPlus_& rhs)
{
	ZPlus_ retval(0);
	// hold enough space 
	retval.digits_.resize(lhs.digits_.size() + rhs.digits_.size(), 0ul);
	// accumulate one-digit multiples
	for (size_t ir = 0; ir < rhs.digits_.size(); ++ir)
		for (size_t il = 0; il < lhs.digits_.size(); ++il)
			retval.digits_[ir + il] += rhs.digits_[ir] * lhs.digits_[il];
	// sweep clean and drop zeroes
	while(retval.Sweep());
	while (!retval.digits_.empty() && !retval.digits_.back())
		retval.digits_.pop_back();
	return retval;
}

ostream& operator<<(ostream& dst, const ZPlus_& n)
{
	for (auto pd = n.digits_.rbegin(); pd != n.digits_.rend(); ++pd)
		dst << *pd;
	return dst;
}

int main(int argc, char* argv[])
{
	int p2 = 1;
	ZPlus_ n(2ul);
	for (int ii = 0; ii < 7; ++ii)
	{
		p2 *= 2;
		n = n * n;
		cout << "2^" << p2 << " = " << n << "\n";
	}
	return 0;
}
2^2 = 4
2^4 = 16
2^8 = 256
2^16 = 65536
2^32 = 4294967296
2^64 = 18446744073709551616
2^128 = 340282366920938463463374607431768211456

Ceylon

"run() is the main function of this module."

shared void run() {

    function multiply(String|Integer|Integer[] top, String|Integer|Integer[] bottom, Integer base = 10) {

        function fromString(String s) =>
                s
                .filter(not(','.equals))
                .map((char) => Integer.parse(char.string))
                .narrow<Integer>()
                .sequence()
                .reversed;

        function toString(Integer[] ints) =>
                ""
                .join(ints.interpose(',', 3))
                .reversed
                .trimLeading((char) => char in "0,");

        function fromInteger(Integer int) => fromString(int.string);

        function convertArg(String|Integer|Integer[] arg) =>
                switch(arg)
                case (is String) fromString(arg)
                case (is Integer) fromInteger(arg)
                case (is Integer[]) arg;

        value a = convertArg(top);
        value b = convertArg(bottom);

        value p = a.size;
        value q = b.size;
        value product = Array.ofSize(p + q, 0);

        for (bIndex->bDigit in b.indexed) {
            variable value carry = 0;
            for (aIndex->aDigit in a.indexed) {
                assert (exists prodDigit = product[aIndex + bIndex]);
                value temp =  prodDigit + carry + aDigit * bDigit;
                carry = temp / base;
                product[aIndex + bIndex] = temp % base;
            }
            assert (exists lastDigit = product[bIndex + p]);
            product[bIndex + p] = lastDigit + carry;
        }

        return toString(product.sequence());
    }

    value twoToThe64th = "18,446,744,073,709,551,616";
    value expectedResult = "340,282,366,920,938,463,463,374,607,431,768,211,456";
    value result = multiply(twoToThe64th, twoToThe64th);

    print("The expected result is ``expectedResult``");
    print("The actual result is   ``result``");
    print("Do they match? ``expectedResult == result then "Yes!" else "No!"``");
}

COBOL

       identification division.
       program-id. long-mul.
       data division.
       replace ==ij-lim== by ==7== ==ir-lim== by ==14==.
       working-storage section.
       1 input-string pic x(26) value "18,446,744,073,709,551,616".
       1 a-table.
        2 a pic 999 occurs ij-lim.
       1 b-table.
        2 b pic 999 occurs ij-lim.
       1 ir-table value all "0".
        2 occurs ij-lim.
         3 ir pic 999 occurs ir-lim.
       1 s-table value all "0".
        2 s pic 999 occurs ir-lim.
       1 display.
        2 temp-result pic 9(6) value 0.
        2 carry pic 999 value 0.
        2 remain pic 999 value 0.
       1 binary.
        2 i pic 9(4) value 0.
        2 j pic 9(4) value 0.
        2 k pic 9(4) value 0.
       procedure division.
       begin.
           move 1 to j
           perform varying i from 1 by 1 until i > ij-lim
               unstring input-string delimited ","
                   into a (i) with pointer j
           end-perform
           move a-table to b-table
           perform intermediate-calc
           perform sum-ir
           perform display-result
       stop run
       .

       intermediate-calc.
           perform varying i from ij-lim by -1 until i < 1
               move 0 to carry
               perform varying j from ij-lim by -1 until j < 1
                   compute temp-result = a (i) * b (j) + carry
                   divide temp-result by 1000 giving carry
                       remainder remain
                   compute k = i + j
                   move remain to ir (i k)
               end-perform
               subtract 1 from k
               move carry to ir (i k)
           end-perform
           .

       sum-ir.
           move 0 to carry
           perform varying k from ir-lim by -1 until k < 1
               move carry to temp-result
               perform varying i from ij-lim by -1 until i < 1
                   compute temp-result = temp-result + ir (i k)
               end-perform
               divide temp-result by 1000 giving carry
                   remainder remain
               move remain to s (k)
           end-perform
           .

       display-result.
           display "   " input-string
           display " * " input-string
           display " = " with no advancing
           perform varying k from 1 by 1
           until k > ir-lim or s (k) not = 0
           end-perform
           if s (k) < 100
               move 1 to i
               inspect s (k) tallying i for leading "0"
               display s (k) (i:) "," with no advancing
               add 1 to k
           end-if
           perform varying k from k by 1 until k > ir-lim
               display s (k) with no advancing
               if k < ir-lim
                   display "," with no advancing
               end-if
           end-perform
           display space
           .

       end program long-mul.
   18,446,744,073,709,551,616
 * 18,446,744,073,709,551,616
 = 340,282,366,920,938,463,463,374,607,431,768,211,456

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

Common Lisp

(defun number->digits (number)
  (do ((digits '())) ((zerop number) digits)
    (multiple-value-bind (quotient remainder) (floor number 10)
      (setf number quotient)
      (push remainder digits))))

(defun digits->number (digits)
  (reduce #'(lambda (n d) (+ (* 10 n) d)) digits :initial-value 0))

(defun long-multiply (a b)
  (labels ((first-digit (list)
             "0 if list is empty, else first element of list."
             (if (endp list) 0
               (first list)))
           (long-add (digitses &optional (carry 0) (sum '()))
             "Do long addition on the list of lists of digits.  Each
              list of digits in digitses should begin with the least
              significant digit.  This is the opposite of the digit
              list returned by number->digits which places the most
              significant digit first.  The digits returned by
              long-add do have the most significant bit first."
             (if (every 'endp digitses)
               (nconc (number->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 (number->digits a)))
         (b (nreverse (number->digits b)))
         (prefix '() (list* 0 prefix))
         (rows '()))
        ((endp b) (digits->number (long-add rows)))
      (let* ((bi (pop b))
             (row (mapcar #'(lambda (ai) (* ai bi)) a)))
        (push (append prefix row) rows)))))
> (long-multiply (expt 2 64) (expt 2 64))
340282366920938463463374607431768211456

Crystal

require "big"

a = 2.to_big_i ** 64

puts "#{a} * #{a} = #{a*a}"
Output:
18446744073709551616 * 18446744073709551616 = 340282366920938463463374607431768211456

D

Using the standard library:

void main() {
    import std.stdio, std.bigint;

    writeln(2.BigInt ^^ 64 * 2.BigInt ^^ 64);
}
Output:
340282366920938463463374607431768211456

Long multiplication, same output:

Translation of: JavaScript
import std.stdio, std.algorithm, std.range, std.ascii, std.string;

auto longMult(in string x1, in string x2) pure nothrow @safe {
    auto digits1 = x1.representation.retro.map!q{a - '0'};
    immutable digits2 = x2.representation.retro.map!q{a - '0'}.array;
    uint[] res;

    foreach (immutable i, immutable d1; digits1.enumerate) {
        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!digits;
    return res.retro.map!(d => digits[d]);
}

void main() {
    immutable two64 = "18446744073709551616";
    longMult(two64, two64).writeln;
}

Dc

This example is incorrect. Please fix the code and remove this message.
Details: Code does not explicitly implement long multiplication

Since Dc has arbitrary precision built-in, the task is no different than a normal multiplication:

2 64^ 2 64^ *p
This example is incorrect. Please fix the code and remove this message.
Details: A Dc solution might be: Represent bignums as numerical strings and implement arithmetic functions on them.

Delphi

Translation of: Go

Copy of core Go answer.

program Long_multiplication;

{$APPTYPE CONSOLE}

uses
  System.SysUtils;

type
  TLongMul = record
  private
    function Add(x, y: TArray<byte>): TArray<byte>;
    function ByteToString(b: TArray<byte>): Ansistring;
    function d(b: byte): Byte;
    function mulDigit(x: TArray<byte>; y: byte): TArray<byte>;
    function mul(x1, y1: AnsiString): AnsiString;
  public
    value: string;
    class operator Multiply(a, b: TLongMul): TLongMul;
    class operator Implicit(a: TLongMul): string;
    class operator Implicit(a: string): TLongMul;
  end;

function TLongMul.d(b: byte): Byte;
begin
  if (b < ord('0')) or (b > ord('9')) then
    raise Exception.Create('digit 0-9 expected: ' + ord(b).ToString);
  Result := ord(b) - ord('0');
end;

class operator TLongMul.Implicit(a: string): TLongMul;
begin
  Result.value := a;
end;

class operator TLongMul.Implicit(a: TLongMul): string;
begin
  Result := a.value;
end;

function TLongMul.Add(x, y: TArray<byte>): TArray<byte>;
begin
  if length(x) < Length(y) then
  begin
    var tmp := y;
    y := x;
    x := tmp;
  end;

  var b: TArray<byte>;
  SetLength(b, length(x) + 1);
  var c: byte := 0;
  for var i := 1 to Length(x) do
  begin
    if i <= Length(y) then
      c := c + d(y[Length(y) - i]);
    var s := d(x[Length(x) - i]) + c;
    c := s div 10;
    b[length(b) - i] := (s mod 10) + ord('0');
  end;
  if c = 0 then
  begin
    Result := b;
    Delete(Result, 0, 1);
    exit;
  end;

  b[0] := c + ord('0');
  Result := b;
end;

function TLongMul.mulDigit(x: TArray<byte>; y: byte): TArray<byte>;
begin
  if y = ord('0') then
  begin
    SetLength(result, 1);
    Result[0] := y;
    exit
  end;

  y := d(y);
  var b: TArray<byte>;
  SetLength(b, length(x) + 1);
  var c: byte := 0;
  for var i := 1 to Length(x) do
  begin
    var s := d(x[Length(x) - i]) * y + c;
    c := s div 10;
    b[length(b) - i] := (s mod 10) + ord('0');
  end;

  if c = 0 then
  begin
    Result := b;
    Delete(Result, 0, 1);
    exit;
  end;

  b[0] := c + ord('0');
  Result := b;
end;

class operator TLongMul.Multiply(a, b: TLongMul): TLongMul;
begin
  Result.value := a.mul(a, b);
end;

function TLongMul.ByteToString(b: TArray<byte>): Ansistring;
begin
  SetLength(Result, length(b));
  move(b[0], Result[1], length(b));
end;

function TLongMul.mul(x1, y1: AnsiString): AnsiString;
var
  x, y: TArray<byte>;
  res: TArray<byte>;
begin
  SetLength(x, length(x1));
  move(x1[1], x[0], length(x1));

  SetLength(y, length(y1));
  move(y1[1], y[0], length(y1));

  res := mulDigit(x, y[length(y) - 1]);

  var zeros: TArray<byte> := [];

  for var i := 2 to Length(y) do
  begin
    SetLength(zeros, Length(zeros) + 1);
    zeros[High(zeros)] := ord('0');

    res := add(res, Concat(mulDigit(x, y[Length(y) - i]), zeros));
  end;

  Result := ByteToString(res);
end;

const
  validate = '340282366920938463463374607431768211456';

var
  num: TLongMul;

begin
  num.value := '18446744073709551616';

  Writeln((num * num).value);
  Writeln(validate);
  Readln;
end.
Output:
340282366920938463463374607431768211456
340282366920938463463374607431768211456

EasyLang

func$ mult a$ b$ .
   a[] = number strchars a$
   b[] = number strchars b$
   len r[] len a[] + len b[]
   for ib = len b[] downto 1
      h = 0
      for ia = len a[] downto 1
         h += r[ia + ib] + b[ib] * a[ia]
         r[ia + ib] = h mod 10
         h = h div 10
      .
      r[ib] += h
   .
   r$ = ""
   for i = 1 to len r[]
      if r$ <> "" or r[i] <> 0 or i = len r[]
         r$ &= r[i]
      .
   .
   return r$
.
print mult "18446744073709551616" "18446744073709551616"

EchoLisp

We implement long multiplication by multiplying polynomials, knowing that the number 1234 is the polynomial x^3 +2x^2 +3x +4 at x=10. As we assume no bigint library is present, long-mul operates on strings.

(lib 'math) ;; for poly multiplication

;; convert string of decimal digits to polynomial
;; "1234" → x^3 +2x^2 +3x +4 
;; least-significant digit first
(define (string->long N)
	(reverse (map string->number (string->list N))))
	
;; convert polynomial to string
(define (long->string N)
(if (pair? N)
   (string-append (number->string (first N)) (long->string (rest N)))  ""))

;; convert poly coefficients to base 10
(define (poly->10 P (carry 0))
(append
	(for/list ((coeff P)) 
		(set! coeff (+ carry coeff )) 
		(set! carry (quotient coeff 10)) ;; new carry
		(modulo coeff 10)) 
	(if(zero? carry) null (list carry)))) ;; remove leading 0 if any

;; long multiplication
;; convert input - strings of decimal digits - to polynomials
;; perform poly multiplication in base 10
;; convert result to string of decimal digits

(define (long-mul A B )
 (long->string (reverse  (poly->10 (poly-mul (string->long A) (string->long B))))))

(define two-64 "18446744073709551616")
(long-mul two-64 two-64)
     "340282366920938463463374607431768211456"

;; check it
(lib 'bigint)
Lib: bigint.lib loaded.
(expt 2 128)
    340282366920938463463374607431768211456

Euphoria

constant base = 1000000000

function atom_to_long(atom a)
    sequence s
    s = {}
    while a>0 do
        s = append(s,remainder(a,base))
        a = floor(a/base)
    end while
    return s
end function

function long_mult(object a, object b)
    sequence c
    if atom(a) then
        a = atom_to_long(a)
    end if
    if atom(b) then
        b = atom_to_long(b)
    end if
    c = repeat(0,length(a)+length(b))
    for i = 1 to length(a) do
        c[i .. i+length(b)-1] += a[i]*b
    end for

    for i = 1 to length(c) do
        if c[i] > base then
            c[i+1] += floor(c[i]/base) -- carry
            c[i] = remainder(c[i],base)
        end if
    end for

    if c[$] = 0 then
        c = c[1..$-1]
    end if
    return c
end function


function long_to_str(sequence a)
    sequence s
    s = sprintf("%d",a[$])
    for i = length(a)-1 to 1 by -1 do
        s &= sprintf("%09d",a[i])
    end for
    return s
end function

sequence a, b, c

a = atom_to_long(power(2,32))
printf(1,"a is %s\n",{long_to_str(a)})

b = long_mult(a,a)
printf(1,"a*a is %s\n",{long_to_str(b)})

c = long_mult(b,b)
printf(1,"a*a*a*a is %s\n",{long_to_str(c)})

Output:

a is 4294967296
a*a is 18446744073709551616
a*a*a*a is 340282366920938463488374607424768211456

F#

This example is incorrect. Please fix the code and remove this message.
Details: The problem is to implement long multiplication, not to demonstrate bignum support.
> let X = 2I ** 64 * 2I ** 64 ;;

val X : System.Numerics.BigInteger = 340282366920938463463374607431768211456

Factor

USING: kernel math sequences ;

: longmult-seq ( xs ys -- zs )
[ * ] cartesian-map
dup length iota [ 0 <repetition> ] map
[ prepend ] 2map
[ ] [ [ 0 suffix ] dip [ + ] 2map ] map-reduce ;

: integer->digits ( x -- xs )  { } swap  [ dup 0 > ] [ 10 /mod swap [ prefix ] dip ] while  drop ;
: digits->integer ( xs -- x )  0 [ swap 10 * + ] reduce ;

: longmult ( x y -- z )  [ integer->digits ] bi@ longmult-seq digits->integer ;
( scratchpad ) 2 64 ^ dup longmult .
340282366920938463463374607431768211456
( scratchpad ) 2 64 ^ dup * .
340282366920938463463374607431768211456

Fortran

Works with: Fortran version 95 and later
module LongMoltiplication
  implicit none

  type longnum
     integer, dimension(:), pointer :: num
  end type longnum

  interface operator (*)
     module procedure longmolt_ll
  end interface

contains
  
  subroutine longmolt_s2l(istring, num)
    character(len=*), intent(in) :: istring
    type(longnum), intent(out) :: num
    
    integer :: i, l

    l = len(istring)

    allocate(num%num(l))

    forall(i=1:l) num%num(l-i+1) = iachar(istring(i:i)) - 48

  end subroutine longmolt_s2l

  ! this one performs the moltiplication
  function longmolt_ll(a, b) result(c)
    type(longnum) :: c
    type(longnum), intent(in) :: a, b
    
    integer, dimension(:,:), allocatable :: t
    integer :: ntlen, i, j

    ntlen = size(a%num) + size(b%num) + 1
    allocate(c%num(ntlen))
    c%num = 0

    allocate(t(size(b%num), ntlen))
    
    t = 0
    forall(i=1:size(b%num), j=1:size(a%num)) t(i, j+i-1) = b%num(i) * a%num(j)

    do j=2, ntlen    
       forall(i=1:size(b%num)) t(i, j) = t(i, j) + t(i, j-1)/10
    end do

    forall(j=1:ntlen) c%num(j) = sum(mod(t(:,j), 10))

    do j=2, ntlen
       c%num(j) = c%num(j) + c%num(j-1)/10
    end do

    c%num = mod(c%num, 10)
    
    deallocate(t)
  end function longmolt_ll


  subroutine longmolt_print(num)
    type(longnum), intent(in) :: num

    integer :: i, j
    
    do j=size(num%num), 2, -1
       if ( num%num(j) /= 0 ) exit
    end do

    do i=j, 1, -1
       write(*,"(I1)", advance="no") num%num(i)
    end do
  end subroutine longmolt_print

end module LongMoltiplication
program Test
  use LongMoltiplication

  type(longnum) :: a, b, r

  call longmolt_s2l("18446744073709551616", a)
  call longmolt_s2l("18446744073709551616", b)

  r = a * b
  call longmolt_print(r)
  write(*,*)

end program Test

FreeBASIC

' version 08-01-2017
' compile with: fbc -s console

Const As UInteger base_ = 1000000000 ' base 1,000,000,000

Function multiply(a1 As String, b1 As String) As String

    Dim As String a = a1, b = b1
    
    Trim(a) : Trim(b) ' remove spaces 
    If Len(a) = 0 Or Len(b) = 0 Then Return "0"

    If Len(a) + Len(b) > 10000 Then 
        Print "number(s) are to big"
        Sleep 5000,1
        Return ""
    End If
    
    If Len(a) < Len(b) Then
        Swap a, b
    End If

    Dim As ULongInt product
    Dim As UInteger carry, i, m, shift
    Dim As UInteger la = Len(a), lb = Len(b)
    Dim As UInteger la9 = la \ 9 + IIf((la Mod 9) = 0, 0, 1)
    Dim As UInteger lb9 = lb \ 9 + IIf((lb Mod 9) = 0, 0, 1)
    Dim As UInteger arr_a(la9), answer((la9 + lb9) + 2)
    Dim As Integer last = la9

    ' make length a, b a multipy of 9
    a = Right((String(9, "0") + a), la9 * 9)
    b = Right((String(9, "0") + b), lb9 * 9)

    For i = 1 To la9
        arr_a(la9 - i +1) = Val(Mid(a, i * 9 -8, 9))
    Next

    Do
        carry = 0
        m = Val(Mid(b, lb9 * 9 -8, 9))
        For i = 1 To la9
            product = CULngInt(arr_a(i)) * m + answer(i + shift) + carry
            carry = product \ base_
            answer(i + shift) = product - carry * base_
        Next
        If carry <> 0 Then
            last = la9 + shift +1
            answer(last) = carry
        End If
        lb9 = lb9 -1
        shift = shift +1
    Loop Until lb9 = 0

    Dim As String tmp = Str(answer(last))
    last = last -1
    While last > 0
        tmp = tmp + Right(String(9,"0") + Str(answer(last)), 9)
        last = last -1
    Wend

    Return tmp

End Function

' ------=< MAIN >=------

Dim As String a = "2", b = "2", answer
Dim As UInteger i = 1, j

For j = 1 To 7
    answer = multiply(a, b)
    a = answer
    b = answer
    i = i + i
    Print using "2 ^ ### = "; i;
    Print answer
Next

Print
Print "-------------------------------------------------"
Print

a = "2" : b = "1" : answer = ""
For j = 1 To 128
    answer = multiply(a, b)
    b = answer
Next
Print "2 ^ 128 = "; answer

' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
Output:
2 ^   2 = 4
2 ^   4 = 16
2 ^   8 = 256
2 ^  16 = 65536
2 ^  32 = 4294967296
2 ^  64 = 18446744073709551616
2 ^ 128 = 340282366920938463463374607431768211456

-------------------------------------------------

2 ^ 128 = 340282366920938463463374607431768211456

Go

// Long multiplication per WP article referenced by task description.
// That is, multiplicand is multiplied by single digits of multiplier
// to form intermediate results.  Intermediate results are accumulated
// for the product.  Used here is the abacus method mentioned by the
// article, of summing intermediate results as they are produced,
// rather than all at once at the end.
//
// Limitations:  Negative numbers not supported, superfluous leading zeros
// not generally removed.
package main

import "fmt"

// argument validation
func d(b byte) byte {
    if b < '0' || b > '9' {
        panic("digit 0-9 expected")
    }
    return b - '0'
}

// add two numbers as strings
func add(x, y string) string {
    if len(y) > len(x) {
        x, y = y, x
    }
    b := make([]byte, len(x)+1)
    var c byte
    for i := 1; i <= len(x); i++ {
        if i <= len(y) {
            c += d(y[len(y)-i])
        }
        s := d(x[len(x)-i]) + c
        c = s / 10
        b[len(b)-i] = (s % 10) + '0'
    }
    if c == 0 {
        return string(b[1:])
    }
    b[0] = c + '0'
    return string(b)
}

// multipy a number by a single digit
func mulDigit(x string, y byte) string {
    if y == '0' {
        return "0"
    }
    y = d(y)
    b := make([]byte, len(x)+1)
    var c byte
    for i := 1; i <= len(x); i++ {
        s := d(x[len(x)-i])*y + c
        c = s / 10
        b[len(b)-i] = (s % 10) + '0'
    }
    if c == 0 {
        return string(b[1:])
    }
    b[0] = c + '0'
    return string(b)
}

// multiply two numbers as strings
func mul(x, y string) string {
    result := mulDigit(x, y[len(y)-1])
    for i, zeros := 2, ""; i <= len(y); i++ {
        zeros += "0"
        result = add(result, mulDigit(x, y[len(y)-i])+zeros)
    }
    return result
}

// requested output
const n = "18446744073709551616"

func main() {
    fmt.Println(mul(n, n))
}

Output:

340282366920938463463374607431768211456

Haskell

import Data.List (transpose, inits)
import Data.Char (digitToInt)

longmult :: Integer -> Integer -> Integer
longmult x y = foldl1 ((+) . (10 *)) (polymul (digits x) (digits y))

polymul :: [Integer] -> [Integer] -> [Integer]
polymul xs ys =
  sum <$>
  transpose
    (zipWith
       (<>)
       (inits $ repeat 0)
       ((\f x -> fmap ((<$> x) . f)) (*) xs ys))

digits :: Integer -> [Integer]
digits = fmap (fromIntegral . digitToInt) . show

main :: IO ()
main = print $ (2 ^ 64) `longmult` (2 ^ 64)
Output:
340282366920938463463374607431768211456

Icon and Unicon

This example is incorrect. Please fix the code and remove this message.
Details: Task is "Implement long multiplication" not "Multiply two numbers using native operators"

Large integers are native to Icon and Unicon. Neither libraries nor special programming is required.

procedure main()
   write(2^64*2^64)
end
Output:
340282366920938463463374607431768211456

J

Solution:

   digits      =: ,.&.":
   polymult    =: +//.@(*/)
   buildDecimal=: 10x&#.

   longmult=: buildDecimal@polymult&digits

Example:

   longmult~ 2x^64
340282366920938463463374607431768211456

Alternatives:
longmult could have been defined concisely:

longmult=: 10x&#.@(+//.@(*/)&(,.&.":))

Or, of course, the task may be accomplished without the verb definitions:

   10x&#.@(+//.@(*/)&(,.&.":))~2x^64
340282366920938463463374607431768211456

Or using the code (+ 10x&*)/@|. instead of #.:

   (+ 10x&*)/@|.@(+//.@(*/)&(,.&.":))~2x^64
340282366920938463463374607431768211456

Or you could use the built-in language support for arbitrary precision multiplication:

   (2x^64)*(2x^64)
340282366920938463463374607431768211456

Explaining the component verbs:

  • digits translates a number to a corresponding list of digits;
   ,.&.": 123
1 2 3
  • polymult (multiplies polynomials): ref. [1]
   1 2 3 (+//.@(*/)) 1 2 3
1 4 10 12 9
  • buildDecimal (translates a list of decimal digits - possibly including "carry" - to the corresponding extended precision number):
   (+ 10x&*)/|. 1 4 10 12 9
15129

or

   10 #. 1 4 10 12 9
15129

Java

Decimal version

This version of the code keeps the data in base ten. By doing this, we can avoid converting the whole number to binary and we can keep things simple, but the runtime will be suboptimal.

public class LongMult {

	private static byte[] stringToDigits(String num) {
		byte[] result = new byte[num.length()];
		for (int i = 0; i < num.length(); i++) {
			char c = num.charAt(i);
			if (c < '0' || c > '9') {
				throw new IllegalArgumentException("Invalid digit " + c
						+ " found at position " + i);
			}
			result[num.length() - 1 - i] = (byte) (c - '0');
		}
		return result;
	}

	public static String longMult(String num1, String num2) {
		byte[] left = stringToDigits(num1);
		byte[] right = stringToDigits(num2);
		byte[] result = new byte[left.length + right.length];
		for (int rightPos = 0; rightPos < right.length; rightPos++) {
			byte rightDigit = right[rightPos];
			byte temp = 0;
			for (int leftPos = 0; leftPos < left.length; leftPos++) {
				temp += result[leftPos + rightPos];
				temp += rightDigit * left[leftPos];
				result[leftPos + rightPos] = (byte) (temp % 10);
				temp /= 10;
			}
			int destPos = rightPos + left.length;
			while (temp != 0) {
				temp += result[destPos] & 0xFFFFFFFFL;
				result[destPos] = (byte) (temp % 10);
				temp /= 10;
				destPos++;
			}
		}
		StringBuilder stringResultBuilder = new StringBuilder(result.length);
		for (int i = result.length - 1; i >= 0; i--) {
			byte digit = result[i];
			if (digit != 0 || stringResultBuilder.length() > 0) {
				stringResultBuilder.append((char) (digit + '0'));
			}
		}
		return stringResultBuilder.toString();
	}

	public static void main(String[] args) {
		System.out.println(longMult("18446744073709551616",
				"18446744073709551616"));
	}
}

Binary version

This version tries to be as efficient as possible, so it converts numbers into binary before doing any calculations. The complexity is higher because of the need to convert to and from base ten, which requires the implementation of some additional arithmetic operations beyond long multiplication itself.

import java.util.Arrays;

public class LongMultBinary {

	/**
	 * A very basic arbitrary-precision integer class. It only handles
	 * non-negative numbers and doesn't implement any arithmetic not necessary
	 * for the task at hand.
	 */
	public static class MyLongNum implements Cloneable {

		/*
		 * The actual bits of the integer, with the least significant place
		 * first. The biggest native integer type of Java is the 64-bit long,
		 * but since we need to be able to store the result of two digits
		 * multiplied, we have to use the second biggest native type, the 32-bit
		 * int. All numeric types are signed in Java, but we don't want to waste
		 * the sign bit, so we need to take extra care while doing arithmetic to
		 * ensure unsigned semantics.
		 */
		private int[] digits;

		/*
		 * The number of digits actually used in the digits array. Since arrays
		 * cannot be resized in Java, we are better off remembering the logical
		 * size ourselves, instead of reallocating and copying every time we need to shrink.
		 */
		private int digitsUsed;

		@Override
		public MyLongNum clone() {
			try {
				MyLongNum clone = (MyLongNum) super.clone();
				clone.digits = clone.digits.clone();
				return clone;
			} catch (CloneNotSupportedException e) {
				throw new Error("Object.clone() threw exception", e);
			}
		}

		private void resize(int newLength) {
			if (digits.length < newLength) {
				digits = Arrays.copyOf(digits, newLength);
			}
		}

		private void adjustDigitsUsed() {
			while (digitsUsed > 0 && digits[digitsUsed - 1] == 0) {
				digitsUsed--;
			}
		}

		/**
		 * "Short" multiplication by one digit. Used to convert strings to long numbers.
		 */
		public void multiply(int multiplier) {
			if (multiplier < 0) {
				throw new IllegalArgumentException(
						"Signed arithmetic isn't supported");
			}
			resize(digitsUsed + 1);
			long temp = 0;
			for (int i = 0; i < digitsUsed; i++) {
				temp += (digits[i] & 0xFFFFFFFFL) * multiplier;
				digits[i] = (int) temp; // store the low 32 bits
				temp >>>= 32;
			}
			digits[digitsUsed] = (int) temp;
			digitsUsed++;
			adjustDigitsUsed();
		}

		/**
		 * "Short" addition (adding a one-digit number). Used to convert strings to long numbers.
		 */
		public void add(int addend) {
			if (addend < 0) {
				throw new IllegalArgumentException(
						"Signed arithmetic isn't supported");
			}
			long temp = addend;
			for (int i = 0; i < digitsUsed && temp != 0; i++) {
				temp += (digits[i] & 0xFFFFFFFFL);
				digits[i] = (int) temp; // store the low 32 bits
				temp >>>= 32;
			}
			if (temp != 0) {
				resize(digitsUsed + 1);
				digits[digitsUsed] = (int) temp;
				digitsUsed++;
			}
		}

		/**
		 * "Short" division (dividing by a one-digit number). Used to convert numbers to strings.
		 * @param divisor The digit to divide by.
		 * @return The remainder of the division.
		 */
		public int divide(int divisor) {
			if (divisor < 0) {
				throw new IllegalArgumentException(
						"Signed arithmetic isn't supported");
			}
			int remainder = 0;
			for (int i = digitsUsed - 1; i >= 0; i--) {
				long twoDigits = (((long) remainder << 32) | (digits[i] & 0xFFFFFFFFL));
				remainder = (int) (twoDigits % divisor);
				digits[i] = (int) (twoDigits / divisor);
			}
			adjustDigitsUsed();
			return remainder;
		}

		public MyLongNum(String value) {
			// each of our 32-bit digits can store at least 9 decimal digit's worth
			this.digits = new int[value.length() / 9 + 1]; 
			this.digitsUsed = 0;
			// To lower the number of bignum operations, handle nine digits at a time.
			for (int i = 0; i < value.length(); i+=9) {
				String chunk = value.substring(i, Math.min(i+9, value.length()));
				int multiplier = 1;
				int addend = 0;
				for (int j=0; j<chunk.length(); j++) {
					char c = chunk.charAt(j);
					if (c < '0' || c > '9') {
						throw new IllegalArgumentException("Invalid digit " + c
								+ " found in input");
					}
					multiplier *= 10;
					addend *= 10;
					addend += c - '0';
				}
				multiply(multiplier);
				add(addend);
			}
		}

		@Override
		public String toString() {
			if (digitsUsed == 0) {
				return "0";
			}
			MyLongNum dummy = this.clone();
			StringBuilder resultBuilder = new StringBuilder(digitsUsed * 9);
			while (dummy.digitsUsed > 0) {
				// To limit the number of bignum divisions, handle nine digits at a time.
				int decimalDigits = dummy.divide(1000000000);
				for (int i=0; i<9; i++) {
					resultBuilder.append((char) (decimalDigits % 10 + '0'));
					decimalDigits /= 10;
				}
			}
			// Trim any leading zeros we may have created.
			while (resultBuilder.charAt(resultBuilder.length()-1) == '0') {
				resultBuilder.deleteCharAt(resultBuilder.length()-1);
			}
			return resultBuilder.reverse().toString();
		}

		/**
		 * Long multiplication.
		 */
		public void multiply(MyLongNum multiplier) {
			MyLongNum left, right;
			// Make sure the shorter number is on the right-hand side to make things a bit more efficient.
			if (this.digitsUsed > multiplier.digitsUsed) {
				left = this;
				right = multiplier;
			} else {
				left = multiplier;
				right = this;
			}
			int[] newDigits = new int[left.digitsUsed + right.digitsUsed];
			for (int rightPos = 0; rightPos < right.digitsUsed; rightPos++) {
				long rightDigit = right.digits[rightPos] & 0xFFFFFFFFL;
				long temp = 0;
				for (int leftPos = 0; leftPos < left.digitsUsed; leftPos++) {
					temp += (newDigits[leftPos + rightPos] & 0xFFFFFFFFL);
					temp += rightDigit * (left.digits[leftPos] & 0xFFFFFFFFL);
					newDigits[leftPos + rightPos] = (int) temp;
					temp >>>= 32;
				}
				// Roll forward any carry we may have.
				int destPos = rightPos + digitsUsed;
				while (temp != 0) {
					temp += (newDigits[destPos] & 0xFFFFFFFFL);
					newDigits[destPos] = (int) temp;
					temp >>>= 32;
					destPos++;
				}
			}
			this.digits = newDigits;
			this.digitsUsed = newDigits.length;
			adjustDigitsUsed();
		}
	}

	public static void main(String[] args) {
		MyLongNum one = new MyLongNum("18446744073709551616");
		MyLongNum two = one.clone();
		one.multiply(two);
		System.out.println(one);
	}

}

JavaScript

Iterative

With integer expression inputs at this scale, JavaScript still gives a slightly lossy result, despite the subsequent digit by digit string concatenation approach.

The problem is that the JavaScript Math.pow expressions become lossy at around 2^54, and Math.pow(2, 64) evaluates to a rounded:

18446744073709552000 rather than the full 18446744073709551616

This means that to handle larger inputs, the multiplication function needs to have string parameters:

function mult(strNum1,strNum2){

    var a1 = strNum1.split("").reverse();
    var a2 = strNum2.toString().split("").reverse();
    var aResult = new Array;
 
    for ( var iterNum1 = 0; iterNum1 < a1.length; iterNum1++ ) {
        for ( var iterNum2 = 0; iterNum2 < a2.length; iterNum2++ ) {
            var 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] %= 10;
            }
        }
    }
    return aResult.reverse().join("");
}


mult('18446744073709551616', '18446744073709551616')
Output:
340282366920938463463374607431768211456

Functional (ES 5)

The function below accepts integer string or native integer arguments, but as JavaScript (unlike Haskell and Python, for example), lacks an arbitrary precision integer type, larger inputs to this function (beyond the scale of c. 2^54) need to take the form of integer strings, to avoid rounding.

For the same reason, the output always takes the form of an arbitrary precision integer string, rather than a native integer data type. (See the largeIntegerString() helper function below)

(function () {
    'use strict';

    // Javascript lacks an unbounded integer type
    // so this multiplication function takes and returns
    // long integer strings rather than any kind of native integer

    // longMult :: (String | Integer) -> (String | Integer) -> String
    function longMult(num1, num2) {
        return largeIntegerString(
            digitProducts(digits(num1), digits(num2))
        );
    }

    // digitProducts :: [Int] -> [Int] -> [Int]
    function digitProducts(xs, ys) {
        return multTable(xs, ys)
            .map(function (zs, i) {
                return Array.apply(null, Array(i))
                    .map(function () {
                        return 0;
                    })
                    .concat(zs);
            })
            .reduce(function (a, x) {
                if (a) {
                    var lng = a.length;

                    return x.map(function (y, i) {
                        return y + (i < lng ? a[i] : 0);
                    })

                } else return x;
            })
    }

    // largeIntegerString :: [Int] -> String
    function largeIntegerString(lstColumnValues) {
        var dctProduct = lstColumnValues
            .reduceRight(function (a, x) {
                var intSum = x + a.carried,
                    intDigit = intSum % 10;

                return {
                    digits: intDigit
                        .toString() + a.digits,
                    carried: (intSum - intDigit) / 10
                };
            }, {
                digits: '',
                carried: 0
            });

        return (dctProduct.carried > 0 ? (
            dctProduct.carried.toString()
        ) : '') + dctProduct.digits;
    }

    // multTables :: [Int] -> [Int] -> [[Int]]
    function multTable(xs, ys) {
        return ys.map(function (y) {
            return xs.map(function (x) {
                return x * y;
            })
        });
    }

    // digits :: (Integer | String) -> [Integer]
    function digits(n) {
        return (typeof n === 'string' ? n : n.toString())
            .split('')
            .map(function (x) {
                return parseInt(x, 10);
            });
    }

    // TEST showing that larged bounded integer inputs give only rounded results
    // whereas integer string inputs allow for full precision on this scale (2^128)

    return {
        fromIntegerStrings: longMult(
            '18446744073709551616',
            '18446744073709551616'
        ),
        fromBoundedIntegers: longMult(
            18446744073709551616,
            18446744073709551616
        )
    };
})();
Output:
{"fromIntegerStrings":"340282366920938463463374607431768211456", 
"fromBoundedIntegers":"340282366920938477630474056040704000000"}

jq

Works with: jq version 1.4

Since the task description mentions 2^64, the following includes "long_power(i)" for computing n^i.

# multiply two decimal strings, which may be signed (+ or -)
def long_multiply(num1; num2):

   def stripsign:
     .[0:1] as $a
     | if $a == "-" then [ -1, .[1:]] 
     elif $a == "+" then [  1, .[1:]] 
     else [1, .]
     end;

  def adjustsign(sign):
     if sign == 1 then . else "-" + . end;

  # mult/2 assumes neither argument has a sign
  def mult(num1;num2):
      (num1 | explode | map(.-48) | reverse) as $a1
    | (num2 | explode | map(.-48) | reverse) as $a2
    | reduce range(0; num1|length) as $i1
        ([];  # result
         reduce range(0; num2|length) as $i2 (.;
  	  ($i1 + $i2) as $ix
  	  | ( $a1[$i1] * $a2[$i2] +
                (if $ix >= length then 0
                 else .[$ix]
                 end) ) as $r
            | if $r > 9 # carrying
              then
                .[$ix + 1] = ($r / 10 | floor) +
                   (if $ix + 1 >= length then 0
                    else .[$ix + 1]
                    end)
                | .[$ix] = $r - ( $r / 10 | floor ) * 10
              else
                .[$ix] = $r
              end
         )
       ) 
    | reverse | map(.+48) | implode;

     (num1|stripsign) as $a1
   | (num2|stripsign) as $a2
   | if $a1[1] == "0" or  $a2[1] == "0" then "0"
     elif $a1[1] == "1" then $a2[1]|adjustsign( $a1[0] * $a2[0] )
     elif $a2[1] == "1" then $a1[1]|adjustsign( $a1[0] * $a2[0] )
     else mult($a1[1]; $a2[1]) | adjustsign( $a1[0] * $a2[0] )
     end;
# Emit (input)^i where input and i are non-negative decimal integers,
# represented as numbers and/or strings.
def long_power(i):
  def power(i):
    tostring as $self
    | (i|tostring) as $i
    | if   $i == "0" then "1"
      elif $i == "1" then $self
      elif $self == "0" then "0"
      else reduce range(1;i) as $_ ( $self; long_multiply(.; $self) )
      end;

  (i|tonumber) as $i
  | if $i < 4 then power($i)
    else ($i|sqrt|floor) as $j
    | ($i - $j*$j) as $k
    | long_multiply( power($j) | power($j) ; power($k) )
  end ;

Example:

 2 | long_power(64) | long_multiply(.;.)
Output:
$ jq -n -f Long_multiplication.jq
"340282366920938463463374607431768211456"

Julia

Works with: Julia version 0.6
Translation of: Python

Module:

module LongMultiplication

using Compat

function addwithcarry!(r, addend, addendpos)
    while true
        pad = max(0, addendpos - lastindex(r))
        append!(r, fill(0, pad))
        addendrst = addend + r[addendpos]
        addend, r[addendpos] = divrem(addendrst, 10)
        iszero(addend) && break
        addendpos += 1
    end
    return r
end

function longmult(mult1::AbstractVector{T}, mult2::AbstractVector{T}) where T <: Integer
    r = T[]
    for (offset1, digit1) in enumerate(mult1), (offset2, digit2) in zip(eachindex(mult2) + offset1 - 1, mult2)
        single_multrst = digits(digit1 * digit2)
        for (addoffset, rstdigit) in zip(eachindex(single_multrst) + offset2 - 1, single_multrst)
            addwithcarry!(r, rstdigit, addoffset)
        end
    end
    return r
end

function longmult(a::T, b::T)::T where T <: Integer
    mult1 = digits(a)
    mult2 = digits(b)
    r = longmult(mult1, mult2)
    return sum(d * T(10) ^ (e - 1) for (e, d) in enumerate(r))
end

function longmult(a::AbstractString, b::AbstractString)
    if !ismatch(r"^\d+", a) || !ismatch(r"^\d+", b)
        throw(ArgumentError("string must contain only digits"))
    end
    mult1 = reverse(collect(Char, a) .- '0')
    mult2 = reverse(collect(Char, b) .- '0')
    r = longmult(mult1, mult2)
    return reverse(join(r))
end

end  # module LongMultiplication

Main:

@show LongMultiplication.longmult(big(2) ^ 64, big(2) ^ 64)
@show LongMultiplication.longmult("18446744073709551616", "18446744073709551616")
Output:
LongMultiplication.longmult(big(2) ^ 64, big(2) ^ 64) = 340282366920938463463374607431768211456
LongMultiplication.longmult("18446744073709551616", "18446744073709551616") = "340282366920938463463374607431768211456"

Kotlin

Translation of: Java
fun String.toDigits() = mapIndexed { i, c ->
    if (!c.isDigit())
        throw IllegalArgumentException("Invalid digit $c found at position $i")
    c - '0'
}.reversed()

operator fun String.times(n: String): String {
    val left = toDigits()
    val right = n.toDigits()
    val result = IntArray(left.size + right.size)

    right.mapIndexed { rightPos, rightDigit ->
        var tmp = 0
        left.indices.forEach { leftPos ->
            tmp += result[leftPos + rightPos] + rightDigit * left[leftPos]
            result[leftPos + rightPos] = tmp % 10
            tmp /= 10
        }
        var destPos = rightPos + left.size
        while (tmp != 0) {
            tmp += (result[destPos].toLong() and 0xFFFFFFFFL).toInt()
            result[destPos] = tmp % 10
            tmp /= 10
            destPos++
        }
    }

    return result.foldRight(StringBuilder(result.size), { digit, sb ->
        if (digit != 0 || sb.length > 0) sb.append('0' + digit)
        sb
    }).toString()
}

fun main(args: Array<out String>) {
    println("18446744073709551616" * "18446744073709551616")
}

Lambdatalk

Natural positive numbers are defined as strings, for instance 123 -> "123". 
{lambda talk} has a small set of primitives working on strings, [equal?, empty?, chars, charAt, substring]

1) helper functions

{def lastchar
 {lambda {:w}
  {charAt {- {chars :w} 1} :w}
}}     
{def butlast
 {lambda {:w}
  {substring 0 {- {chars :w} 1} :w}
}} 
{def zeros
 {lambda {:n}
  {if {< :n 1}
   then 
   else 0{zeros {- :n 1}}
}}}

2) add function

{def add 
 {def add.r
 {lambda {:a :b :c :d}
  {if {equal? :a #}
   then {if {equal? :d 1} then 1 else}{butlast :c}
   else {let { {:a :a} {:b :b} {:c :c}
               {:d {+ :d {lastchar :a} {lastchar :b} }} }
    {add.r {butlast :a} {butlast :b} {lastchar :d}:c
         {if {equal? {chars :d} 1} then 0 else 1}}
 }}}}
 {lambda {:a :b}
  {{lambda {:a :b :n}
    {add.r #{zeros {- :n {chars :a}}}:a
           #{zeros {- :n {chars :b}}}:b  # 0}
  } :a :b {max {chars :a} {chars :b}}}
}}  

3) mul function

{def mul
 {def muln
  {lambda {:a :b :n}
   {if {< :n 1}
    then :b
    else {muln :a {add :a :b} {- :n 1}}
 }}}
 {def mul.r
  {lambda {:a :b :c :n}
   {if {equal? :b #}
    then :c
    else {mul.r :a {butlast :b}
         {add {muln :a 0 {lastchar :b}}{zeros :n} :c} {+ :n 1}}
 }}}
 {lambda {:a :b}
  {mul.r :a #:b 0 0}
}} 

4) applying to the task

Due to JS numbers limits, we compute first 2^32 using the JS pow function, then 2^64 and 2^128 using the mul function.

2^32 = '{def p32 {pow 2 32}}          -> '{p32}  = 4294967296
2^64 = '{def p64 {mul {p32} {p32}}}   -> '{p64}  = 18446744073709551616
2^128 = '{def p128 {mul {p64} {p64}}} -> '{p128} = 340282366920938463463374607431768211456 

5) a more effective implementation 

Lambdatalk can be helped by the lib_BN javascript library from Jonas Raoni Soares Silva
and stored in a wiki page called by a {require lib_BN} command, computing becomes fast:  

2^32 = {def p32 {BN.pow 2 32}}        -> {p32}  = 4294967296
2^64 = {def p64 {BN.* {p32} {p32}}}   -> {p64}  = 18446744073709551616
2^128 = {def p128 {BN.* {p64} {p64}}} -> {p128} = 340282366920938463463374607431768211456

This can be tested in http://lambdaway.free.fr/lambdaspeech/?view=numbers8

Liberty BASIC

Works with: Just BASIC
Works with: Run BASIC
'[RC] long multiplication

'now, count 2^64
print "2^64"
a$="1"
for i = 1 to 64
    a$ = multByD$(a$, 2)
next
print a$
print "(check with native LB)"
print 2^64
print "(looks OK)"

'now let's do b$*a$ stuff
print
print "2^64*2^64"
print longMult$(a$, a$)
print "(check with native LB)"
print 2^64*2^64
print "(looks OK)"

end
'---------------------------------------
function longMult$(a$, b$)
    signA = 1
    if left$(a$,1) = "-" then a$ = mid$(a$,2): signA = -1
    signB = 1
    if left$(b$,1) = "-" then b$ = mid$(b$,2): signB = -1

    c$ = ""
    t$ = ""
    shift$ = ""
    for i = len(a$) to 1 step -1
        d = val(mid$(a$,i,1))
        t$ = multByD$(b$, d)
        c$ = addLong$(c$, t$+shift$)
        shift$ = shift$ +"0"
    'print d, t$, c$ 
    next
    if signA*signB<0 then c$ = "-" + c$
    'print c$
    longMult$ = c$
end function

function multByD$(a$, d)
'multiply a$ by digit d
c$ = ""
carry = 0
for i = len(a$) to 1 step -1
        a = val(mid$(a$,i,1))
        c = a*d+carry
        carry = int(c/10)
        c = c mod 10
        'print a, c
        c$ = str$(c)+c$ 
next
    if carry>0 then c$ = str$(carry)+c$
    'print c$
    multByD$ = c$
end function

function addLong$(a$, b$)
'add a$ + b$, for now only positive
    l = max(len(a$), len(b$))
    a$=pad$(a$,l)
    b$=pad$(b$,l)
    c$ = "" 'result
    carry = 0
    for i = l to 1 step -1
        a = val(mid$(a$,i,1))
        b = val(mid$(b$,i,1))
        c = a+b+carry
        carry = int(c/10)
        c = c mod 10
        'print a, b, c
        c$ = str$(c)+c$
    next
    if carry>0 then c$ = str$(carry)+c$
    'print c$
    addLong$ = c$
end function

function pad$(a$,n)  'pad from right with 0 to length n
     pad$ = a$
     while len(pad$)<n
        pad$ = "0"+pad$
     wend
end function

Lobster

Translation of: Java
Translation of Java binary version, but with base 1000000000
import std

// Very basic arbitrary-precision integers
// - only non-negative numbers
// - doesn't implement any arithmetic not necessary for the task at hand...

let base = 1000000000

class Bign:
    digits: [int] // little endian, of base base
    digitsUsed: int

    def clone():
        return Bign { digits: copy(digits), digitsUsed: digitsUsed }

    def resize(newLength):
        while digits.length < newLength:
            digits.push(0)

    def adjustDigitsUsed():
        while digitsUsed > 0 and digits[digitsUsed - 1] == 0:
                digitsUsed -= 1

    // multiplication by one digit; used to convert string to Bign
    def muldigit(multiplier : int):
        if (multiplier < 0):
            return // "Signed arithmetic isn't supported"
        resize(digitsUsed + 1)
        var temp = 0
        for(digitsUsed) i:
            temp += digits[i] * multiplier
            digits[i] = temp % base
            temp /= base
        digits[digitsUsed] = temp
        digitsUsed += 1
        adjustDigitsUsed()

    // addition of one digit; used to convert string to Bign
    def adddigit(addend: int):
        if (addend < 0):
            return // "Signed arithmetic isn't supported"
        var temp = addend
        var i = 0
        while i < digitsUsed and temp != 0:
            temp += digits[i]
            digits[i] = temp % base
            temp /= base
            i += 1
        if temp != 0:
            resize(digitsUsed + 1)
            digits[digitsUsed] = temp
            digitsUsed += 1

    def bign2str():
        var i = digitsUsed
        if i == 0:
            return "0"
        i -= 1
        var s = string(digits[i])
        while i > 0:
            i -= 1
            s += number_to_string(digits[i], 10, 9)
        return s

def str2bign(value):
    // each of our Bign digits can store 9 decimal digits
    let this = Bign { digits: map(value.length() / 9 + 1): 0, digitsUsed: 0 }
    // handle nine digits at a time
    var i = 0
    while i < value.length:
        var multiplier = 1
        var addend = 0
        for(min(9, value.length() - i)) j:
            let c = value[i+j]
            //if (c < '0' or c > '9') -- what!?
            multiplier *= 10
            addend *= 10
            addend += c - '0'
        this.muldigit(multiplier)
        this.adddigit(addend)
        i += 9
    return this

// Long multiplication

def bign_multiply(this, multiplier):
    // Make sure the shorter number is on the right side to make things a bit more efficient
    let left  = if (this.digitsUsed > multiplier.digitsUsed): this else: multiplier
    let right = if (this.digitsUsed > multiplier.digitsUsed): multiplier else: this
    let newDigits = map(left.digitsUsed + right.digitsUsed): 0
    for(right.digitsUsed) rightPos:
        let rightDigit = right.digits[rightPos]
        var temp = 0
        for(left.digitsUsed) leftPos:
            temp += newDigits[leftPos + rightPos]
            temp += rightDigit * left.digits[leftPos]
            newDigits[leftPos + rightPos] = temp % base
            temp /= base
        // Roll forward any carry we may have
        let destPos = rightPos + left.digitsUsed
        while temp != 0:
            temp += newDigits[destPos]
            newDigits[destPos] = temp % base
            temp /= base
            destPos +- 1
    let bign = Bign { digits: newDigits, digitsUsed: newDigits.length }
    bign.adjustDigitsUsed()
    return bign

let one = str2bign("18446744073709551616")
let two = one.clone()
var pro = one.bign_multiply(two)
print(bign2str(pro))
Output:
340282366920938463463374607431768211456

Maple

longmult := proc(a::integer,b::integer)
    local A,B,m,n,i,j;
    # Note, return a*b; works in Maple for any sized integer
    A := convert(a,base,10);
    B := convert(b,base,10);
    m := numelems(A);
    n := numelems(B);
    add( add( A[i]*B[j]*10^(j-1), j=1..n )*10^(i-1), i=1..m );
end;

> longmult( 2^64, 2^64 );
                    340282366920938463463374607431768211456

Mathematica/Wolfram Language

We define the long multiplication function:

 LongMultiplication[a_,b_]:=Module[{d1,d2},
  d1=IntegerDigits[a]//Reverse;
  d2=IntegerDigits[b]//Reverse;
  Sum[d1[[i]]d2[[j]]*10^(i+j-2),{i,1,Length[d1]},{j,1,Length[d2]}]
 ]

Example:

 n1 = 2^64;
 n2 = 2^64;
 LongMultiplication[n1, n2]

gives back:

 340282366920938463463374607431768211456

To check the speed difference between built-in multiplication (which is already arbitrary precision) we multiply two big numbers (2^8000 has 2409 digits!) and divide their timings:

 n1=2^8000;
 n2=2^8000;
 Timing[LongMultiplication[n1,n2]][[1]]
 Timing[n1 n2][[1]]
 Floor[%%/%]

gives back:

 72.9686
 7.*10^-6
 10424088

So our custom function takes about 73 second, the built-in function a couple of millionths of a second, so the long multiplication is about 10.5 million times slower! Mathematica uses Karatsuba multiplication for large integers, which is several magnitudes faster for really big numbers. Making it able to multiply in about a second; the final result has 9542426 digits; result omitted for obvious reasons.

NetRexx

Translation of: REXX

A reworking of the example at Rexx Version 2.

/* NetRexx */
options replace format comments java crossref symbols nobinary

numeric digits 100

runSample(arg)
return

-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
method multiply(multiplier, multiplicand) public static
  result = ''
  mpa = s2a(multiplier)
  mpb = s2a(multiplicand)
  r_ = 0
  rim = 1
  loop bi = 1 to mpb[0]
    loop ai = 1 to mpa[0]
      ri = ai + bi -1
      p_ = mpa[ai] * mpb[bi]
      loop i_ = ri by 1 until p_ = 0
        s_ = r_[i_] + p_
        r_[i_] = s_ // 10
        p_ = s_ % 10
        end i_
      rim = rim.max(i_)
      end ai
    end bi
  r_[0] = rim
  result = a2s(r_)
  result = result.strip('l', 0)
  if result = '' then result = 0
  return result

-- .............................................................................
-- copy characters of a numeric string into a corresponding array
-- digits are numbered 1 to n from right to left
method s2a(numbr) private static
  result = 0
  lstr = numbr.length()
  loop z_ = 1 to lstr
    result[z_] = numbr.substr(lstr - z_ + 1, 1)
    end z_
  result[0] = lstr
  return result

-- .............................................................................
-- turn the array of digits into a numeric string
method a2s(numbr) private static
  result = ''
  loop z_ = numbr[0] to 1 by -1
    result = result || numbr[z_]
    end z_
  return result

-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
method runSample(arg) private static
  mms = [ -
                      123',  '123, -
                      012',  '12, -
             123456789012' , '44444444444, -
                  2 ** 64' , '2**64, -
                        0'                      ,0 ' -
  ]
  ok = 0
  errors = 0

  loop mm over mms
    parse mm multiplier . ',' multiplicand .
    builtIn = multiplier * multiplicand
    calculated = multiply(multiplier, multiplicand)
    say 'Calculate' multiplier + 0 'x' multiplicand + 0
    say 'Built in:' builtIn
    say 'Derived: ' calculated
    say
    if builtIn = calculated then ok = ok + 1
    else                         errors = errors + 1
    end mm
  say ok 'ok'
  say errors 'not ok'

  return
Output:
Calculate 123 x 123
Built in: 15129
Derived:  15129

Calculate 12 x 12
Built in: 144
Derived:  144

Calculate 123456789012 x 44444444444
Built in: 5486968400478463649328
Derived:  5486968400478463649328

Calculate 18446744073709551616 x 18446744073709551616
Built in: 340282366920938463463374607431768211456
Derived:  340282366920938463463374607431768211456

Calculate 0 x 0
Built in: 0
Derived:  0

5 ok
0 not ok

Nim

Translation of: C
import strutils

proc ti(a: char): int = ord(a) - ord('0')

proc longmulti(a, b: string): string =
  var
    i, j = 0
    k = false

  # either is zero, return "0"
  if a == "0" or b == "0":
    return "0"

  # see if either a or b is negative
  if a[0] == '-':
    i = 1; k = not k
  if b[0] == '-':
    j = 1; k = not k

  # if yes, prepend minus sign if needed and skip the sign
  if i > 0 or j > 0:
    result = if k: "-" else: ""
    result.add longmulti(a[i..a.high], b[j..b.high])
    return

  result = repeat('0', a.len + b.len)

  for i in countdown(a.high, 0):
    var carry = 0
    var k = i + b.len
    for j in countdown(b.high, 0):
      let n = ti(a[i]) * ti(b[j]) + ti(result[k]) + carry
      carry = n div 10
      result[k] = chr(n mod 10 + ord('0'))
      dec k
    result[k] = chr(ord(result[k]) + carry)

  if result[0] == '0':
    result[0..result.high-1] = result[1..result.high]

echo longmulti("-18446744073709551616", "-18446744073709551616")

Output:

3402823669209384634633746074317682114566

Oforth

Oforth handles arbitrary precision integers, so there is no need to implement long multiplication :

Output:
2 64 pow dup * println
340282366920938463463374607431768211456

But, if long multiplication was to be implemented :

A natural is implemented as a list of integers with base 1000000000 (in order to print them easier)

Just multiplication is implemented here.

Number Class new: Natural(v)
 
Natural method: initialize  := v ;
Natural method: _v  @v ;
 
Natural classMethod: newValues super new ;
Natural classMethod: newFrom   asList self newValues ;
 
Natural method: *(n)
| v i j l x k |
   n _v ->v
   ListBuffer initValue(@v size v size + 1+, 0) ->l
 
   v size loop: i [
      i v at dup ->x 0 ifEq: [ continue ]
      0 @v size loop: j [
         i j + 1- ->k
         j @v at x * + l at(k) + 1000000000 /mod k rot l put
         ]
      k 1+ swap l put
      ]
   while(l last 0 == l size 0 <> and) [ l removeLast drop ]
   l dup freeze Natural newValues ;
 
Natural method: <<
| i |
   @v last <<
   @v size 1 - loop: i [ @v at(@v size i -) <<wjp(0, JUSTIFY_RIGHT, 8) ] ;
Output:
>Natural newFrom(2 16 pow) .s
[1] (Natural) 65536
ok
>dup * .s
[1] (Natural) 4294967296
ok
>dup * .s
[1] (Natural) 18446744073709551616
ok
>dup * .s
[1] (Natural) 340282366920938463463374607431768211456
ok
>_v .s
[1] (List) [768211456, 374607431, 938463463, 282366920, 340]
ok

Ol

Ol already supports long numbers "out-of-the-box".

(define x (* 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2)) ; 2^64

(print (* x x))
340282366920938463463374607431768211456

PARI/GP

long(a,b)={
  a=eval(Vec(a));
  b=eval(Vec(b));
  my(c=vector(#a+#b),carry=0);
  for(i=1,#a,
    for(j=1,#b,
      c[i+j]+=a[i]*b[j]
    )
  );
  forstep(i=#c,1,-1,
    c[i] += carry;
    carry = c[i] \ 10;
    c[i] = c[i] % 10
  );
  for(i=1,#c,
    if(c[i], return(concat(apply(s->Str(s),vector(#c+1-i,j,c[i+j-1])))))
  );
  "0"
};
long("18446744073709551616","18446744073709551616")

Output:

%1 = "340282366920938463463374607431768211456"

Pascal

Extracted from a programme to calculate and factor the number (two versions) in Frederick Pohl's book The Gold at the Starbow's End, and compute Godel encodings of text. Compiles with the Free Pascal Compiler. The original would compile with Turbo Pascal (and used pointers to allow access to the "heap" storage scheme) except that does not allow functions to return a "big number" data aggregate, and it is so much nicer to be able to write X:=BigMult(A,B); The original has a special "square" calculation but this task is to exhibit long multiplication. However, raising to a power by iteration is painful, so a special routine for that.

Program TwoUp; Uses DOS, crt;
{Concocted by R.N.McLean (whom God preserve), Victoria university, NZ.}
 Procedure Croak(gasp: string);
  Begin
   Writeln;
   Write(Gasp);
   HALT;
  End;

 const BigBase = 10;		{The base of big arithmetic.}
 const BigEnuff = 333;		{The most storage possible is 65532 bytes with Turbo Pascal.}
 type  BigNumberIndexer = word;	{To access 0:BigEnuff BigNumberDigit data.}
 type  BigNumberDigit = byte;	{The data.}
 type  BigNumberDigit2 = word;	{Capable of digit*digit + carry. Like, 255*255 = 65025}

 type BigNumber =		{All sorts of arrangements are possible.}
  Record				{Could include a sign indication.}
   TopDigit: BigNumberDigit;		{Finger the high-order digit.}
   digit: array[0..BigEnuff] of byte;	{The digits: note the "downto" in BigShow.}
  end;					{Could add fractional digits too. Endless, endless.}

 Procedure BigShow(var a: BigNumber);	{Print the number.}
  var i: integer;	{A stepper.}
  Begin
   for i:=a.TopDigit downto 0 do	{Thus high-order to low, as is the custom.}
    if BigBase = 10 then write(a.digit[i])	{Constant following by the Turbo Pascal compiler}
     else if BigBase = 100 then Write(a.digit[i] div 10,a.digit[i] mod 10)	{Means that there will be no tests.}
      else write(a.digit[i],',');		{And dead code will be omitted.}
  End;

 Procedure BigZero(var A: BigNumber); {A:=0;}
  Begin;
   A.TopDigit:=0;
   A.Digit[0]:=0;
  End;
 Procedure BigOne(var A: BigNumber);  {A:=1;}
  Begin;
   A.TopDigit:=0;
   A.Digit[0]:=1;
  End;
 Function BigInt(n: longint): BigNumber; {A:=N;}
  var l: BigNumberIndexer;
  Begin
   l:=0;
   if n < 0 then croak('Negative integers are not yet considered.');
   repeat		{At least one digit is to be placed.}
    if l > BigEnuff then Croak('BigInt overflowed!');	{Oh dear.}
    BigInt.Digit[l]:=N mod BigBase;	{The low-order digit.}
    n:=n div BigBase;			{Shift down a digit.}
    l:=l + 1;				{Count in anticipation.}
   until N = 0;			{Still some number left?}
   BigInt.TopDigit:=l - 1;	{Went one too far.}
  End;

 Function BigMult(a,b: BigNumber): BigNumber;	{x:=BigMult(a,b);}
{Suppose the digits of A are a5,a4,a3,a2,a1,a0...
 To multiply A and B.
                               a5   a4   a3   a2   a1   a0: six digits, d1
                                x   b4   b3   b2   b1   b0: five digits, d2
                               ---------------------------
                             a5b0 a4b0 a3b0 a2b0 a1b0 a0b0
                        a5b1 a4b1 a3b1 a2b1 a1b1 a0b1
                   a5b2 a4b2 a3b2 a2b2 a1b2 a0b2
              a5b3 a4b3 a3b3 a2b3 a1b3 a0b3
         a5b4 a4b4 a3b4 a2b4 a1b4 a0b4
   -------------------------------------------------------
   carry    9    8    7    6    5    4    3    2    1    0: at least nine digits,
   -------------------------------------------------------  = d1 + d2 - 1
   But the indices are also the powers, so the highest power is 9 = 5 + 4,
and a possible tenth for any carry.}
  var X: BigNumber;		{Scratchpad, so b:=BigMult(a,b); doesn't overwrite b as it goes...}
  var d: BigNumberDigit;	{A digit.}
  var c: BigNumberDigit;	{A carry.}
  var dd: BigNumberDigit2;	{A digit product.}
  var i,j,l: BigNumberIndexer;	{Steppers.}
  Begin
   if ((A.TopDigit = 0) and (A.Digit[0] = 0))
    or((B.TopDigit = 0) and (B.Digit[0] = 0)) then begin BigZero(BigMult); exit; end;
   l:=A.TopDigit + B.TopDigit;       {Minimal digit requirement. (Counting is from zero)}
   if l > BigEnuff then Croak('BigMult will overflow.');
   for i:=l downto 0 do X.Digit[i]:=0;	{Clear for action.}
   for i:=0 to A.TopDigit do		{Arbitrarily, choose A on the one hand.}
    begin				{Though there could be a better choice.}
     d:=A.Digit[i];			{Select the digit.}
     if d <> 0 then			{What the hell. One in BigBase chance.}
      begin				{But not this time.}
       l:=i;				{Locate the power of BigBase.}
       c:=0;				{Start this digit's multiply pass.}
       for j:=0 to B.TopDigit do	{Stepping along B's digits.}
        begin				{One by one.}
         dd:=BigNumberDigit2(B.Digit[j])*d + X.Digit[l] + c;	{The deed.}
         X.Digit[l]:=dd mod BigBase;	{Place the new digit.}
         c:=dd div BigBase;		{And extract the carry.}
         l:=l + 1;			{Ready for the next power up.}
        end;				{Advance to it.}
       if c > 0 then			{The multiply done, place the carry.}
        begin				{Ah. We *will* use the next power up.}
         if l > BigEnuff then Croak('BigMultX has overflowed.');	{Oh dear.}
         X.Digit[l]:=c;		{Thus as if BigMult..Digit[l] was zeroed.}
         l:=l + 1;			{Preserve the one-too-far for the last case}
        end;				{So much for a carry at the end of a pass.}
      end;				{So much for a non-zero digit.}
    end;			{On to another digit to multiply with.}
   X.TopDigit:=l - 1;	{Remember the one-too-far.}
   BigMult:=X;		{Deliver, possibly scragging A or B, or, both!}
 End; {of BigMult.}

 Procedure BigPower(var X: BigNumber; P: longint); {Replaces X by X**P}
  var A,W: BigNumber;	{Scratchpads}
  label up;
  Begin		{Each squaring doubles the power, melding nicely with binary reduction.}
   if P <= 0 then Croak('Negative powers are not accommodated!');
   BigOne(A);		{x**0 = 1}
   W:=X;		{Holds X**1, 2, 4, 8, etc.}
up:if P mod 2 = 1 then A:=BigMult(A,W);	{Bit on, so include this order.}
   P:=P div 2;		{Halve the power contrariwise to W's doubling.}
   if P > 0 then 	{Still some power to come?}
    begin		{Yes.}
     W:=BigMult(W,W);	{Step up to the next bit's power.}
     goto up;		{And see if it is "on".}
    end;		{Odd layout avoids multiply testing P > 0.}
   X:=A;		{The result.}
  End;

 var X: BigNumber;
 var p: longint;
 BEGIN
  ClrScr;
  WriteLn('To calculate  x = 2**64, then x*x via multi-digit long multiplication.');
  p:=64;		{As per the specification.}
  X:=BigInt(2);		{Start with 2.}
  BigPower(X,p);	{First stage: 2**64}
  Write ('x = 2**',p,' = '); BigShow(X);
  WriteLn;
  X:=BigMult(X,X);	{Second stage.}
  Write ('x*x = ');BigShow(X);	{Can't have Write('x*x = ',BigShow(BigMult(X,X))), after all. Oh well.}
 END.

Output:

To calculate  x = 2**64, then x*x via multi-digit long multiplication.
x = 2**64 = 18446744073709551616
x*x = 340282366920938463463374607431768211456

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

Phix

base 10^9

Translation of: Euphoria

Simple longhand multiplication. To keep things as simple as possible, this does not handle negative numbers.
If bcd1 is a number split into digits 0..9, bcd9 is a number split into "digits" 000,000,000..999,999,999, which fit in an integer.
They are held lsb-style mainly so that trimming a trailing 0 does not alter their value.

constant base = 1_000_000_000
 
function bcd9_mult(sequence a, sequence b)
 
    sequence c = repeat(0,length(a)+length(b))
    for i=1 to length(a) do
        integer j = i+length(b)-1
        c[i..j] = sq_add(c[i..j],sq_mul(a[i],b))
    end for
 
    for i=1 to length(c) do
        atom ci = c[i]
        if ci>base then
            c[i+1] += floor(ci/base) -- carry
            c[i] = remainder(ci,base)
        end if
    end for
 
    if c[$]=0 then
        c = c[1..$-1]
    end if
    return c
end function
 
function atom_to_bcd9(atom a)
sequence s = {}
    while a>0 do
        s = append(s,remainder(a,base))
        a = floor(a/base)
    end while
    return s
end function
 
function bcd9_to_str(sequence a)
string s = sprintf("%d",a[$])
    for i=length(a)-1 to 1 by -1 do
        s &= sprintf("%09d",a[i])
    end for
    -- (might want to trim leading 0s here)
    return s
end function
 
sequence a, b, c
 
a = atom_to_bcd9(power(2,32))
printf(1,"a is %s\n",{bcd9_to_str(a)})
 
b = bcd9_mult(a,a)
printf(1,"a*a is %s\n",{bcd9_to_str(b)})
 
c = bcd9_mult(b,b)
printf(1,"a*a*a*a is %s\n",{bcd9_to_str(c)})
Output:
a is 4294967296
a*a is 18446744073709551616
a*a*a*a is 340282366920938463488374607488768211456

string

with javascript_semantics
function mul(string a, b)
    bool bSign = false
    if a[1]='-' then {bSign,a} = {not bSign, a[2..$]} end if
    if b[1]='-' then {bSign,b} = {not bSign, b[2..$]} end if
    string res = repeat('0',length(a)+length(b))
    --
    -- Note that i,j,k are used as negative indexes, working 
    -- from the right hand least significant digit leftwards.
    --
    for i=1 to length(a) do
        integer j=1, k=i, c=0
        while j<=length(b) or c do
            c += res[-k]-'0'
            if j<=length(b) then
                c += (a[-i]-'0')*(b[-j]-'0')
                j += 1
            end if          
            res[-k] = remainder(c,10)+'0'
            c = floor(c/10)
            k += 1
        end while
    end for
    res = trim_head(res,"0")
    if bSign then res = '-'&res end if
    return res
end function
?mul("18446744073709551616","18446744073709551616")
Output:
"340282366920938463488374607488768211456"

builtin

Library: Phix/mpfr

(same output as immediately above)

include mpfr.e
mpz a = mpz_init("18446744073709551616") -- or:
--mpz a = mpz_init(); mpz_ui_pow_ui(a,2,64)
mpz_mul(a,a,a)
?mpz_get_str(a)

PHP

<?php
function longMult($a, $b)
{
  $as = (string) $a;
  $bs = (string) $b;
  for($pi = 0, $ai = strlen($as) - 1; $ai >= 0; $pi++, $ai--)
    {
      for($p = 0; $p < $pi; $p++)
        {
          $regi[$ai][] = 0;
        }
      for($bi = strlen($bs) - 1; $bi >= 0; $bi--)
        {
          $regi[$ai][] = $as[$ai] * $bs[$bi];
        }
    }
  return $regi;
}

function longAdd($arr)
{
  $outer = count($arr);
  $inner = count($arr[$outer-1]) + $outer;
  for($i = 0; $i <= $inner; $i++)
    {
      for($o = 0; $o < $outer; $o++)
        {
          $val  = isset($arr[$o][$i]) ? $arr[$o][$i] : 0;
          @$sum[$i] += $val;
        }
    }
  return $sum;
}

function carry($arr)
{
  for($i = 0; $i < count($arr); $i++)
    {
      $s = (string) $arr[$i];
      switch(strlen($s))
        {
          case 2:
            $arr[$i] = $s{1};
            @$arr[$i+1] += $s{0};
            break;
          case 3:
            $arr[$i] = $s{2};
            @$arr[$i+1] += $s{0}.$s{1};
            break;
        }
    }
  return ltrim(implode('',array_reverse($arr)),'0');
}

function lm($a,$b)
{
  return carry(longAdd(longMult($a,$b)));
}

if(lm('18446744073709551616','18446744073709551616') == '340282366920938463463374607431768211456')
  {
    echo 'pass!';
  }; // 2^64 * 2^64

PicoLisp

(de multi (A B)
   (setq A (format A) B (reverse (chop B)))
   (let Result 0
      (for (I . X) B 
         (setq Result (+ Result (* (format X) A (** 10 (dec I)))))) ) )

PL/I

/* Multiply a by b, giving c. */
multiply: procedure (a, b, c);
   declare (a, b, c) (*) fixed decimal (1);
   declare (d, e, f) (hbound(a,1)) fixed decimal (1);
   declare pr (-hbound(a,1) : hbound(a,1)) fixed decimal (1);
   declare p fixed decimal (2), (carry, s) fixed decimal (1);
   declare neg bit (1) aligned;
   declare (i, j, n, offset) fixed binary (31);

   n = hbound(a,1);
   d = a;
   e = b;
   s = a(1) + b(1);
   neg = (s = 9);
   if a(1) = 9 then call complement (d);
   if b(1) = 9 then call complement (e);
   pr = 0;
   offset = 0; carry = 0;
   do i = n to 1 by -1;
      do j = n to 1 by -1;
         p = d(i) * e(j) + pr(j-offset) + carry;
         if p > 9 then do; carry = p/10; p = mod(p, 10); end; else carry = 0;
         pr(j-offset) = p;
      end;
      offset = offset + 1;
   end;
   do i = hbound(a,1) to 1 by -1;
      c(i) = pr(i);
   end;
   do i = -hbound(a,1) to 1;
      if pr(i) ^= 0 then signal fixedoverflow;
   end;
   if neg then call complement (c);
end multiply;

complement: procedure (a);
   declare a(*) fixed decimal (1);
   declare i fixed binary (31), carry fixed decimal (1);
   declare s fixed decimal (2);

   carry = 1;
   do i = hbound(a,1) to 1 by -1;
      s = 9 - a(i) + carry;
      if s > 9 then do; s = s - 10; carry = 1; end; else carry = 0;
      a(i) = s;
   end;
end complement;

Calling sequence:

   a = 0; b = 0; c = 0;
   a(60) = 1;
   do i = 1 to 64; /* Generate 2**64 */
      call add (a, a, b);
      put skip;
      call output (b);
      a = b;
   end;
   call multiply (a, b, c);
   put skip;
   call output (c);

Final output:

18446744073709551616
 340282366920938463463374607431768211456

PL/M

Based on the Algol W sample, Uses bytes instead of integers to hold the digits. Ony handles positive numbers.

100H: /* LONG MULTIPLICATION OF LARGE INTEGERS                               */
   /* LARGE INTEGERS ARE REPRESENTED BY ARRAYS OF BYTES WHOSE VALUES ARE     */
   /* A SINGLE DECIMAL DIGIT OF THE NUMBER                                   */
   /* THE LEAST SIGNIFICANT DIGIT OF THE LARGE INTEGER IS IN ELEMENT 1       */
   /* ELEMENT 0 CONTAINS THE NUMBER OF DIGITS THE NUMBER HAS                 */
   BDOS: PROCEDURE( FN, ARG ); /* CP/M BDOS SYSTEM CALL                      */
      DECLARE FN BYTE, ARG ADDRESS;
      GOTO 5;
   END BDOS;
   PRINT$CHAR:   PROCEDURE( C ); DECLARE C BYTE;    CALL BDOS( 2, C ); END;
   PRINT$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
   DECLARE PRINT$NL LITERALLY 'PRINT$STRING( .( 0DH, 0AH, ''$'' ) )';

   DECLARE LONG$INTEGER  LITERALLY '(201)BYTE';
   DECLARE DIGIT$BASE    LITERALLY '10';

   /* PRINTS A LONG INTEGER                                                  */
   PRINT$LONG$INTEGER: PROCEDURE( N$PTR );
      DECLARE N$PTR ADDRESS;
      DECLARE N BASED N$PTR LONG$INTEGER;
      DECLARE ( D, F ) BYTE;
      F = N( 0 );
      DO D = 1 TO N( 0 );
         CALL PRINT$CHAR( N( F ) + '0' );
         F = F - 1;
      END;
   END PRINT$LONG$INTEGER;
   /* IMPLEMENTS LONG MULTIPLICATION, C IS SET TO A * B                      */
   /*     C CAN BE THE SAME LONG$INTEGER AS A OR B                           */
   LONG$MULTIPLY: PROCEDURE( A$PTR, B$PTR, C$PTR );
      DECLARE ( A$PTR, B$PTR, C$PTR ) ADDRESS;
      DECLARE ( A BASED A$PTR, B BASED B$PTR, C BASED C$PTR ) LONG$INTEGER;
      DECLARE MRESULT LONG$INTEGER;
      DECLARE RPOS    BYTE;

      /* MULTIPLIES THE LONG INTEGER IN B BY THE INTEGER A, THE RESULT       */
      /*     IS ADDED TO C, STARTING FROM DIGIT START                        */
      /*     OVERFLOW IS IGNORED                                             */
      MULTIPLY$ELEMENT: PROCEDURE( A, B$PTR, C$PTR, START );
         DECLARE ( B$PTR, C$PTR )                 ADDRESS;
         DECLARE ( A, START )                     BYTE;
         DECLARE ( B BASED B$PTR, C BASED C$PTR ) LONG$INTEGER;
         DECLARE ( CDIGIT, D$CARRY, BPOS, CPOS )  BYTE;
         D$CARRY = 0;
         CPOS    = START;
         DO BPOS = 1 TO B( 0 );
            CDIGIT = C( CPOS ) + ( A * B( BPOS ) ) + D$CARRY;
            IF CDIGIT < DIGIT$BASE THEN D$CARRY = 0;
            ELSE DO;
               /* HAVE DIGITS TO CARRY                                       */
               D$CARRY = CDIGIT  /  DIGIT$BASE;
               CDIGIT  = CDIGIT MOD DIGIT$BASE;
            END;
            C( CPOS ) = CDIGIT;
            CPOS = CPOS + 1;
         END;
         C( CPOS ) = D$CARRY;
         /* REMOVE LEADING ZEROS BUT IF THE NUMBER IS 0, KEEP THE FINAL 0    */
         DO WHILE( CPOS > 1 AND C( CPOS ) = 0 );
            CPOS = CPOS - 1;
         END;
         C( 0 ) = CPOS;
      END MULTIPLY$ELEMENT ;

      /* THE RESULT WILL BE COMPUTED IN MRESULT, ALLOWING A OR B TO BE C     */
      DO RPOS = 1 TO LAST( MRESULT ); MRESULT( RPOS ) = 0; END;
      /* MULTIPLY BY EACH DIGIT AND ADD TO THE RESULT                        */
      DO RPOS = 1 TO A( 0 );
         IF A( RPOS ) <> 0 THEN DO;
            CALL MULTIPLY$ELEMENT( A( RPOS ), B$PTR, .MRESULT, RPOS );
         END;
      END;
      /* RETURN THE RESULT IN C                                              */
      DO RPOS = 0 TO MRESULT( 0 ); C( RPOS ) = MRESULT( RPOS ); END;
   END;

   /* CALCULATE AND OUTPUT 2^128                                             */
   DECLARE ( TWO$TO$64, TWO$TO$128 ) LONG$INTEGER;
   DECLARE ( PWR, TPOS ) BYTE;
   /* CONSTRUCT 2^64 IN TWO$TO$64                                            */
   DO TPOS = 0 TO LAST( TWO$TO$64 ); TWO$TO$64( TPOS ) = 0; END;
   TWO$TO$64( 0 ) = 1;
   TWO$TO$64( 1 ) = 2;
   PWR          = 1;
   DO WHILE PWR < 64;
      CALL LONG$MULTIPLY( .TWO$TO$64, .TWO$TO$64, .TWO$TO$64 );
      PWR  = PWR + PWR;
   END;
   /* CONSTRUCT 2^128                                                        */
   TWO$TO$128( 0 ) = 1;
   TWO$TO$128( 1 ) = 0;
   CALL LONG$MULTIPLY( .TWO$TO$64, .TWO$TO$64, .TWO$TO$128 );
   CALL PRINT$STRING( .( '2', 05EH, '128: $' ) ); /* 05EH = "^" IN ASCII     */
   CALL PRINT$LONG$INTEGER( .TWO$TO$128 );
   CALL PRINT$NL;
EOF
Output:
2^128: 340282366920938463463374607431768211456

PowerShell

Implementation

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

Library Method

Works with: PowerShell version 4.0
[BigInt]$n = [Math]::Pow(2,64)
[BigInt]::Multiply($n,$n)

Output:

 
340282366920938463463374607431768211456

Prolog

Arbitrary precision arithmetic is native in most Prolog implementations.

 ?- X is 2**64 * 2**64.
X = 340282366920938463463374607431768211456.

PureBasic

Explicit Implementation

Structure decDigitFmt ;decimal digit format
  Array Digit.b(0) ;contains each digit of number, right-most digit is index 0
  digitCount.i ;zero based
  sign.i ; {x < 0} = -1, {x = 0} = 0,  {x > 0} = 1
EndStructure

Global zero_decDigitFmt.decDigitFmt ;represents zero in the decimal digit format

;converts string representation of integer into the digit format, number can include signus but no imbedded spaces
Procedure stringToDecDigitFmt(numString.s, *x.decDigitFmt)
  Protected *c.Character, digitIdx, digitCount
  If numString And *x
    *c.Character = @numString
    Repeat
      Select *c\c
        Case '0' To '9', '-', '+'
          *c + SizeOf(Character)
        Default
          numString = Left(numString, *c - @numString)
          Break
      EndSelect
    ForEver
    *c = @numString
    Select  *c\c
      Case '-'
        *x\sign = -1
        *c + SizeOf(Character)
      Case '+'
        *x\sign = 1
        *c + SizeOf(Character)
      Case '0' To '9'
        *x\sign = 1
    EndSelect
    
    numString = LTrim(PeekS(*c), "0") ;remove leading zeroes
    If numString = "" ;is true if equal to zero or if only a signus is present
      CopyStructure(@zero_decDigitFmt, *x, decDigitFmt)
      ProcedureReturn
    EndIf
    *c = @numString
    
    digitCount = Len(PeekS(*c)) - 1
    Dim *x\Digit(digitCount)
    *x\digitCount = digitCount
    
    digitIdx = 0
    While *c\c
      If *c\c >= '0' And *c\c <= '9'
        *x\Digit(digitCount - digitIdx) = *c\c - '0'
        digitIdx + 1
        *c + SizeOf(Character)
      Else
        Break
      EndIf
    Wend 
    
  EndIf 
EndProcedure

;converts digit format representation of integer into string representation
Procedure.s decDigitFmtToString(*x.decDigitFmt)
  Protected i, number.s
  If *x
    If *x\sign = 0
      number = "0"
    Else
        
      For i = *x\digitCount To 0 Step -1
        number + Str(*x\Digit(i))
      Next
      number = LTrim(number, "0")
      If *x\sign = -1
        number = "-" + number
      EndIf 
    EndIf 
  EndIf
  
  ProcedureReturn number
EndProcedure

;handles only positive numbers and zero, negative numbers left as an exercise for the reader ;)
Procedure add_decDigitFmt(*a.decDigitFmt, *b.decDigitFmt, *sum.decDigitFmt, digitPos = 0) ;*sum contains the result of (*a ) * 10^digitPos + (*b)
  Protected carry, i, newDigitCount, workingSum, a_dup.decDigitFmt
  
  If *a And *b And *sum
    
    If *a = *sum: CopyStructure(*a, @a_dup, decDigitFmt): *a = @a_dup: EndIf ;handle special case of  *sum + *b = *sum
    If *b <> *sum: CopyStructure(*b, *sum, decDigitFmt): EndIf ;handle general case of *a + *b = *sum and special case of *a + *sum = *sum
    
    ;calculate number of digits needed for sum and resize array of digits if necessary
    newDigitCount = *a\digitCount + digitPos
    If newDigitCount >= *sum\digitCount
      If *sum\digitCount = newDigitCount And *sum\Digit(*sum\digitCount) <> 0
        newDigitCount + 1
      EndIf 
      
      If *sum\digitCount <> newDigitCount
        *sum\digitCount = newDigitCount
        Redim *sum\Digit(*sum\digitCount)
      EndIf
    EndIf
    
    i = 0 
    Repeat
      If i <= *a\digitCount
        workingSum = *a\Digit(i) + *sum\Digit(digitPos) + carry
      Else
        workingSum = *sum\Digit(digitPos) + carry
      EndIf 
      
      If workingSum > 9
        carry = 1
        workingSum - 10
      Else
        carry = 0
      EndIf
      *sum\Digit(digitPos)  = workingSum
      digitPos + 1
      i + 1
    Until i > *a\digitCount And carry = 0
    
    If *a\sign <> 0 Or *sum\sign <> 0
      *sum\sign = 1 ;only handle positive numbers and zero for now
    EndIf 
  EndIf
EndProcedure

Procedure multiply_decDigitFmt(*a.decDigitFmt, *b.decDigitFmt, *product.decDigitFmt) ;*product contains the result of (*a) * (*b)
  Protected i, digitPos, productSignus
  Protected Dim multTable.decDigitFmt(9)
  Protected NewList digitProduct.decDigitFmt()
  
  If *a And *b And *product
    If *a\sign = 0 Or *b\sign = 0
      CopyStructure(zero_decDigitFmt, *product, decDigitFmt)
      ProcedureReturn
    EndIf 
    
    If *b\digitCount > *a\digitCount: Swap *a, *b: EndIf 
    
    ;build multiplication table
    CopyStructure(*a, @multTable(1), decDigitFmt): multTable(1)\sign = 1 ;always positive
    For i = 2 To 9
      add_decDigitFmt(*a, multTable(i - 1), multTable(i))
    Next 
    
    ;collect individual digit products for later summation; these could also be added as we go along
    For i = 0 To *b\digitCount
      AddElement(digitProduct())
      digitProduct() = multTable(*b\Digit(i))
    Next
    
    ;determine sign of product
    If *a\sign <> *b\sign
      productSignus = -1
    Else
      productSignus = 1
    EndIf 
    
    digitPos = 0
    CopyStructure(zero_decDigitFmt, *product, decDigitFmt)
    ForEach digitProduct()
      add_decDigitFmt(digitProduct(), *product, *product, digitPos)
      digitPos + 1
    Next
    *product\sign = productSignus ;set sign of product
  EndIf
EndProcedure

;handles only positive integer exponents or an exponent of zero, does not raise an error for 0^0
Procedure exponent_decDigitFmt(*a.decDigitFmt, exponent, *product.decDigitFmt)
  Protected i, a_dup.decDigitFmt
  If *a And *product And exponent >= 0
    If *a = *product: CopyStructure(*a, @a_dup, decDigitFmt): *a = @a_dup: EndIf
    stringToDecDigitFmt("1", *product)
    For i = 1 To exponent: multiply_decDigitFmt(*product, *a, *product): Next
  EndIf
EndProcedure

If OpenConsole()
  Define a.decDigitFmt, product.decDigitFmt
  
  stringToDecDigitFmt("2", a)
  exponent_decDigitFmt(a, 64, a) ;2^64
  multiply_decDigitFmt(a, a, product)
  PrintN("The result of 2^64 * 2^64 is " + decDigitFmtToString(product))
  Print(#crlf$ + #crlf$ + "Press ENTER to exit"): Input()
  CloseConsole()
EndIf

Output:

The result of 2^64 * 2^64 is 340282366920938463463374607431768211456

Library Method

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.

XIncludeFile "decimal.pbi"

Define.Decimal *a, *b
*a=PowerDecimal(IntegerToDecimal(2),IntegerToDecimal(64))
*b=TimesDecimal(*a,*a,#NoDecimal)

Print("2^64*2^64 = "+DecimalToString(*b))

Outputs

2^64*2^64 = 340282366920938463463374607431768211456

Python

(Note that Python comes with arbitrary length integers).

#!/usr/bin/env python
print 2**64*2**64
Works with: Python version 3.0
Translation of: Perl
#!/usr/bin/env python

def add_with_carry(result, addend, addendpos):
    while True:
        while len(result) < addendpos + 1:
            result.append(0)
        addend_result = str(int(addend) + int(result[addendpos]))
        addend_digits = list(addend_result)
        result[addendpos] = addend_digits.pop()

        if not addend_digits:
            break
        addend = addend_digits.pop()
        addendpos += 1

def longhand_multiplication(multiplicand, multiplier):
    result = []
    for multiplicand_offset, multiplicand_digit in enumerate(reversed(multiplicand)):
        for multiplier_offset, multiplier_digit in enumerate(reversed(multiplier), start=multiplicand_offset):
            multiplication_result = str(int(multiplicand_digit) * int(multiplier_digit))

            for addend_offset, result_digit_addend in enumerate(reversed(multiplication_result), start=multiplier_offset):
                add_with_carry(result, result_digit_addend, addend_offset)

    result.reverse()

    return ''.join(result)

if __name__ == "__main__":
    sixtyfour = "18446744073709551616"

    onetwentyeight = longhand_multiplication(sixtyfour, sixtyfour)
    print(onetwentyeight)

Shorter version:

Translation of: Haskell
Works with: Python version 3.7
'''Long multiplication'''

from functools import reduce


def longmult(x, y):
    '''Long multiplication.'''
    return reduce(
        digitSum,
        polymul(digits(x), digits(y)), 0
    )


def digitSum(a, x):
    '''Left to right decimal digit summing.'''
    return a * 10 + x


def polymul(xs, ys):
    '''List of specific products.'''
    return map(
        lambda *vs: sum(filter(None, vs)),
        *[
            [0] * i + zs for i, zs in
            enumerate(mult_table(xs, ys))
        ]
    )


def mult_table(xs, ys):
    '''Rows of all products.'''
    return [[x * y for x in xs] for y in ys]


def digits(x):
    '''Digits of x as a list of integers.'''
    return [int(c) for c in str(x)]


if __name__ == '__main__':
    print(
        longmult(2 ** 64, 2 ** 64)
    )


Quackery

Long multiplication as it was taught at primary school, using natural numbers (including zero) in the base 10 numeral system, with a slight variation in that it maintains a running total rather than adding up the intermediate results column-wise at the end. Starting point is "learn your tables". Presumptions are an understanding of "equal to" and "not equal to", and the ability to split a one or two digit number into tens and units.

(Splitting a number into tens and units is achieved using /mod to keep the tables compact. Numbers in the addition and multiplication tables could be represented as, for example, [ 8 1 ] instead of 81 and 10 /mod replaced with unpack or, as the nests only contains numbers, do.)

In addition to the specified task, we were always encouraged to show our workings.

( ------------- preamble to task, some i/o related words ------------- )

  [ [] swap witheach
      [ char 0 -
        swap join ] ]                    is $->long    (     $ --> L   )

  [ number$ $->long ]                    is long       (     n --> L   )
  
  [ reverse behead 
    swap witheach
      [ swap 10 * + ] ]                  is long->num  (     L --> n   )

  [ reverse
    witheach echo ]                      is echolong   (     L -->     )

( ------------------------- task starts here ------------------------- )

  [ [ table
      [ 0  1  2  3  4  5  6  7  8  9 ]
      [ 1  2  3  4  5  6  7  8  9 10 ]
      [ 2  3  4  5  6  7  8  9 10 11 ]
      [ 3  4  5  6  7  8  9 10 11 12 ]
      [ 4  5  6  7  8  9 10 11 12 13 ]
      [ 5  6  7  8  9 10 11 12 13 14 ]
      [ 6  7  8  9 10 11 12 13 14 15 ]
      [ 7  8  9 10 11 12 13 14 15 16 ]
      [ 8  9 10 11 12 13 14 15 16 17 ]
      [ 9 10 11 12 13 14 15 16 17 18 ] ]
    swap peek 10 /mod ]                  is add        (   n n --> n n )

  [ dip add add dip [ add nip ] swap ]   is addc       ( n n c --> n c )

  [ over size
    over size -
    dup dip
      [ 0 < if swap ]
    abs times
      [ 0 join ] ]                       is zeropad    (   L L --> L L )

  [ zeropad            ( when adding two numbers of different lengths ) 
    0 temp put         ( leading zeroes are added to make the lengths )
    [] unrot witheach  ( equal. This is implicit when the calculation )
      [ dip behead     ( done by hand, and performed by zeropad here. )
        temp take
        addc
        temp put
        rot swap join swap ]
    drop
    temp take dup 0 !=
    iff join else drop ]                 is longadd    (   L L --> L   )

  [ [ table
      [ 0  0  0  0  0  0  0  0  0  0 ]
      [ 0  1  2  3  4  5  6  7  8  9 ]
      [ 0  2  4  6  8 10 12 14 16 18 ]
      [ 0  3  6  9 12 15 18 21 24 27 ]
      [ 0  4  8 12 16 20 24 28 32 36 ]
      [ 0  5 10 15 20 25 30 35 40 45 ]
      [ 0  6 12 18 24 30 36 42 48 54 ]
      [ 0  7 14 21 28 35 42 49 56 63 ]
      [ 0  8 16 24 32 40 48 56 64 72 ]
      [ 0  9 18 27 36 45 54 63 72 81 ] ]
    swap peek 10 /mod ]                  is mult       (   n n --> n n )

  [ dip mult add dip [ add nip ] swap ]  is multc      ( n n c --> n c )

  [ dup 0 = iff 
      [ 2drop 0 long ] done
    0 temp put
    [] unrot swap witheach
      [ over temp take
        multc
        temp put
        swap dip join ]
    drop
    temp take dup 0 !=
    iff join else drop ]                 is shortmult  (   L n --> L   )

  [ dup 0 long != iff
    [ 0 swap join ] ]                    is timesten   (     L --> L   )

  [ dup 0 long = iff 
      [ 2drop 0 long ] done
    0 long unrot
    witheach
      [ dip dup shortmult
        rot longadd swap
        timesten ] 
     drop ]                              is longmult   (   L L --> L   )

( ------------------------ additional to task ------------------------ )

  [ stack ]                              is linelength (       --> s   )

  [ linelength share times
      [ char - emit ]
    cr ]                                 is separator  (       -->     )

  [ linelength share
    over size - times sp
    echolong cr ]                        is showlong   (     L -->     )

  [ over size
    over size + linelength put
    over showlong
    dup showlong
    separator
    dup 0 long = iff 
      [ 2drop 0 long ] done
    0 long unrot
    witheach
      [ dip dup shortmult
        dup showlong
        rot longadd swap
        timesten ] 
    drop
    separator
    showlong
    separator
    linelength release ]                 is workings   (   L L -->     )

( --------------------------- demonstration -------------------------- )

say "Using long multiplication: "
2  64 ** long dup longmult dup echolong cr

say "Using built-in arithmetic: "
2 128 ** dup echo cr cr

swap long->num = iff
  say "10/10, Gold star!"
else
  say "0/10, See me after class."

cr cr
say "(Show your workings.)" cr cr
2 64 ** long dup workings cr
Output:
Using long multiplication: 340282366920938463463374607431768211456
Using built-in arithmetic: 340282366920938463463374607431768211456

10/10, Gold star!

(Show your workings.)

                    18446744073709551616
                    18446744073709551616
----------------------------------------
                   110680464442257309696
                   184467440737095516160
                 11068046444225730969600
                 18446744073709551616000
                922337203685477580800000
               9223372036854775808000000
             166020696663385964544000000
                                       0
           12912720851596686131200000000
           55340232221128654848000000000
         1291272085159668613120000000000
                                       0
        73786976294838206464000000000000
       737869762948382064640000000000000
     12912720851596686131200000000000000
    110680464442257309696000000000000000
    737869762948382064640000000000000000
   7378697629483820646400000000000000000
 147573952589676412928000000000000000000
 184467440737095516160000000000000000000
----------------------------------------
 340282366920938463463374607431768211456
----------------------------------------

R

Using GMP

Library: gmp
library(gmp)
a <- as.bigz("18446744073709551616")
mul.bigz(a,a)
"340282366920938463463374607431768211456"

A native implementation

This code is more verbose than necessary, for ease of understanding.

longmult <- function(xstr, ystr)
{
   #get the number described in each string
   getnumeric <- function(xstr) as.numeric(unlist(strsplit(xstr, "")))
   
   x <- getnumeric(xstr)
   y <- getnumeric(ystr)
   
   #multiply each pair of digits together
   mat <- apply(x %o% y, 1, as.character)
   
   #loop over columns, then rows, adding zeroes to end of each number in the matrix to get the correct positioning
   ncols <- ncol(mat)
   cols <- seq_len(ncols)
   for(j in cols)
   {
      zeroes <- paste(rep("0", ncols-j), collapse="") 
      mat[,j] <- paste(mat[,j], zeroes, sep="")  
   }
   
   nrows <- nrow(mat)
   rows <- seq_len(nrows)
   for(i in rows)
   {
      zeroes <- paste(rep("0", nrows-i), collapse="") 
      mat[i,] <- paste(mat[i,], zeroes, sep="")  
   }
   
   #add zeroes to the start of the each number, so they are all the same length
   len <- max(nchar(mat))
   strcolumns <- formatC(cbind(as.vector(mat)), width=len)
   strcolumns <- gsub(" ", "0", strcolumns)
   
   #line up all the numbers below each other
   strmat <- matrix(unlist(strsplit(strcolumns, "")), byrow=TRUE, ncol=len)
   
   #convert to numeric and add them
   mat2 <- apply(strmat, 2, as.numeric)
   sum1 <- colSums(mat2)
   
   #repeat the process on each of the totals, until each total is a single digit
   repeat
   {
      ntotals <- length(sum1)
      totals <- seq_len(ntotals)
      for(i in totals)
      {
         zeroes <- paste(rep("0", ntotals-i), collapse="")
         sum1[i] <- paste(sum1[i], zeroes, sep="")
      }
      len2 <- max(nchar(sum1))
      strcolumns2 <- formatC(cbind(as.vector(sum1)), width=len2)
      strcolumns2 <- gsub(" ", "0", strcolumns2)
      strmat2 <- matrix(unlist(strsplit(strcolumns2, "")), byrow=TRUE, ncol=len2)
      mat3 <- apply(strmat2, 2, as.numeric)
      sum1 <- colSums(mat3)
      if(all(sum1 < 10)) break
   }
   
   #Concatenate the digits together
   ans <- paste(sum1, collapse="")
   ans
}

a <- "18446744073709551616"
longmult(a, a)
"340282366920938463463374607431768211456"

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

Raku

(formerly Perl 6)

Works with: rakudo version 2015-09-17

For efficiency (and novelty), this program explicitly implements long multiplication, but in base 10000. That base was chosen because multiplying two 5-digit numbers can overflow a 32-bit integer, but two 4-digit numbers cannot.

sub num_to_groups ( $num ) { $num.flip.comb(/.**1..4/)».flip     };
sub groups_to_num ( @g   ) { [~] flat @g.pop, @g.reverse».fmt('%04d') };

sub long_multiply ( Str $x, Str $y ) {
    my @group_sums;
    for flat 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;
}

my $str = '18446744073709551616';
long_multiply( $str, $str ).say;

# cross-check with native implementation
say +$str * +$str;
Output:
340282366920938463463374607431768211456
340282366920938463463374607431768211456

REXX

version 1

This REXX version supports:

  •   leading signs
  •   decimal points
  •   automatically adjusting the number of decimal digits needed

Programming note:   &&   is REXX's   exclusive or   operand.

/*REXX program  performs  long multiplication  on  two numbers  (without the "E").      */
numeric digits 300                               /*be able to handle gihugeic input #s. */
parse arg x y .                                  /*obtain optional arguments from the CL*/
if x=='' | x==","  then x= 2**64                 /*Not specified?  Then use the default.*/
if y=='' | y==","  then y= x                     /* "      "         "   "   "     "    */
if x<0  &&  y<0    then sign= '-'                /*there only a single negative number? */
                   else sign=                    /*no, then result sign must be positive*/
xx=x;    x=strip(x, 'T', .);      x1= left(x, 1) /*remove any trailing decimal points.  */
yy=y;    y=strip(y, 'T', .);      y1= left(y, 1) /*   "    "     "        "       "     */
if x1=='-' | x1=="+"  then x= substr(x, 2)       /*remove a leading  ±  sign.           */
if y1=='-' | y1=="+"  then y= substr(y, 2)       /*   "   "    "     "    "             */
parse var x  '.' xf;  parse var y  "." yf        /*obtain the fractional part of X and Y*/
#= length(xf || yf)                              /*#: digits past the decimal points (.)*/
x= space( translate( x, , .),  0)                /*remove decimal point if there is any.*/
y= space( translate( y, , .),  0)                /*   "       "     "    "    "   "  "  */
Lx= length(x);  Ly=length(y)                     /*get the lengths of the new  X and Y. */
numeric digits max(digits(), Lx + Ly)            /*use a new  decimal digits  precision.*/
$= 0                                             /*$:  is the product  (so far).        */
                  do j=Ly  by -1  for Ly         /*almost like REXX does it, ··· but no.*/
                  $= $  +  ((x*substr(y, j, 1))copies(0, Ly-j) )
                  end   /*j*/
f= length($) - #                                 /*does product has enough decimal digs?*/
if f<0  then $=copies(0, abs(f) + 1)$            /*Negative?  Add leading 0s for INSERT.*/
say 'long mult:'  xx  "*"  yy  '──►'   sign || strip( insert(., $, length($) - #), 'T', .)
say ' built─in:'  xx  "*"  yy  '──►'   xx*yy     /*stick a fork in it,  we're all done. */
output   when using the default inputs:
long mult: 18446744073709551616 * 18446744073709551616 ──► 340282366920938463463374607431768211456
 built─in: 18446744073709551616 * 18446744073709551616 ──► 340282366920938463463374607431768211456
output   when using the input of:     123   -456789000
long mult: 123 * -456789000 ──► -56185047000
 built─in: 123 * -456789000 ──► -56185047000
output   when using the input of:     -123.678   +456789000
long mult: -123.678 * +456789000 ──► -56494749942.000
 built─in: -123.678 * +456789000 ──► -56494749942.000

version 2

/* REXX **************************************************************
* While REXX can multiply arbitrary large integers 
* here is the algorithm asked for by the task description
* 13.05.2013 Walter Pachl
*********************************************************************/
cnt.=0
Numeric Digits 100
Call test 123 123
Call test 12 12
Call test 123456789012 44444444444
Call test 2**64 2**64
Call test 0 0
say cnt.0ok 'ok'
say cnt.0nok 'not ok'
Exit
test:
  Parse Arg a b
  soll=a*b
  haben=multiply(a b)
  Say 'soll =' soll
  Say 'haben=' haben
  If haben<>soll Then 
    cnt.0nok=cnt.0nok+1
  Else
    cnt.0ok=cnt.0ok+1
  Return

multiply: Procedure
/* REXX **************************************************************
* Multiply(a b) -> a*b
*********************************************************************/
  Parse Arg a b
  Call s2a 'a'
  Call s2a 'b'
  r.=0
  rim=1
  r0=0
  Do bi=1 To b.0
    Do ai=1 To a.0
      ri=ai+bi-1
      p=a.ai*b.bi
      Do i=ri by 1 Until p=0
        s=r.i+p
        r.i=s//10
        p=s%10
        End
      rim=max(rim,i)
      End
    End
  res=strip(a2s('r'),'L','0')
  If res='' Then
    res='0'
  Return res

s2a:
/**********************************************************************
* copy characters of a string into a corresponding array
* digits are numbered 1 to n fron right to left
**********************************************************************/
  Parse arg name
  string=value(name)
  lstring=length(string)
  do z=1 to lstring
    Call value name'.'z,substr(string,lstring-z+1,1)
    End
  Call value name'.0',lstring
  Return

a2s:
/**********************************************************************
* turn the array of digits into a string
**********************************************************************/
  call trace 'o'
  Parse Arg name
  ol=''
  Do z=rim To 1 By -1
    ol=ol||value(name'.z')
    End
  Return ol

Output:

soll = 15129
haben= 15129
soll = 144
haben= 144
soll = 5486968400478463649328
haben= 5486968400478463649328
soll = 340282366920938463463374607431768211456
haben= 340282366920938463463374607431768211456
soll = 0
haben= 0
5 ok
0 not ok

Ring

This example is incorrect. Please fix the code and remove this message.
Details: Task is "Implement long multiplication" not "Multiply two numbers using native operators"
decimals(0)
see pow(2,64)*pow(2,64) + nl

Output:

340282366920938463463374607431768211456

RPL

To solve the task, we need to develop part of a BigInt library to handle very long integers as strings. Addition has been optimized with a 10-digit batch size, but multiplication is long (and slow), as required.

RPL code Comment
≪ → a 
  ≪ 1 WHILE a OVER DUP SUB "0" == REPEAT 1 + END
     a SWAP OVER SIZE SUB
≫ ≫ ‘NoZero’ STO

≪ WHILE DUP REPEAT "0" ROT + SWAP 1 - END DROP
≫ ‘ZeroFill' STO

≪ DUP2 SIZE SWAP SIZE → sb sa
  ≪ IF sa sb < THEN SWAP END 
     sa sb - ABS 
≫ ≫ 'Swap→Zero' STO

≪ Swap→Zero ZeroFill → a b
  ≪ "" 1 CF 
     a SIZE 1 FOR j 
       1 FS?C 
       a j 9 - j SUB STR→ + b j 9 - j SUB STR→ + 
       IF DUP 1E10 ≥ THEN 1E10 - 1 SF END 
       →STR 10 OVER SIZE - ZeroFill SWAP + 
    -10 STEP 
    1 FS? →STR SWAP + NoZero 
≫ ≫ ‘ADDbig’ STO

≪ → a d
  ≪ "" 0 
     a SIZE 1 FOR j
       a j DUP SUB STR→ d * + 
       10 MOD LAST / IP 
       SWAP →STR ROT + SWAP
    -1 STEP 
    IF THEN LAST →STR SWAP + END
≫ ≫ ‘DigitMul’ STO

≪ Swap→Zero DROP → a b
  ≪ "0" b SIZE 1 FOR j 
        a b j DUP SUB STR→ DigitMul 
        "" b SIZE j - ZeroFill + ADDbig
     -1 STEP
≫ ≫ ‘MULbig’ STO
NoZero ( "0..0xyz" → "xyz" ) 
count leading zeros 
keep the rest


ZeroFill ( "xyz" n → "0..0xyz" ) 

 
Swap→Zero ( a b → a b Δlength ) 
swap a and b if length(a) < length(b)
return length difference


ADDbig ( "a" "b" → "a+b" ) 
res = "" ; carry = 0
for j = length(a) downto 1 step 10
  digits = carry
  digits += a[j-9..j] + b[j-9..j]
  if b > 1E10 then digits -= 1E10 ; carry = 1
  convert digits to 10-char string with leading zeros
next j
prepend carry and remove leading zeros


DigitMul ( "a" d → "a*d" ) 
res = "" ; carry = 0
for j = length(a) downto 1 
   digit = a[j]*d
   carry = digit // 10
   digit %= 10
next j
prepend carry


MULbig ( "a" "b" → "a*b" ) 
sum = "0" ; for j = length(b) downto 1
  tmp = a * b[j]
  shift left tmp length(b)-j times, then add to sum
next j

 "18446744073709551616" DUP MULbig
Output:
1: "340282366920938463463374607431768211456"

Ruby

Translation of: Tcl
def longmult(x,y)
  result = [0]
  j = 0
  y.digits.each do |m|
    c = 0
    i = j
    x.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
 
n=2**64
printf "         %d * %d = %d\n", n, n, n*n
printf "longmult(%d, %d) = %d\n", n, n, longmult(n,n)
         18446744073709551616 * 18446744073709551616 = 340282366920938463463374607431768211456
longmult(18446744073709551616, 18446744073709551616) = 340282366920938463463374607431768211456

Run BASIC

Works with: Just BASIC
Works with: Liberty BASIC
' Same code and result as in Liberty BASIC
Output:
Same as iberty BASIC entry.

Scala

This implementation does not rely on an arbitrary precision numeric type. Instead, only single digits are ever multiplied or added, and all partial results are kept as string.

def addNums(x: String, y: String) = {
  val padSize = x.length max y.length
  val paddedX = "0" * (padSize - x.length) + x
  val paddedY = "0" * (padSize - y.length) + y
  val (sum, carry) = (paddedX zip paddedY).foldRight(("", 0)) {
    case ((dx, dy), (acc, carry)) =>
      val sum = dx.asDigit + dy.asDigit + carry
      ((sum % 10).toString + acc, sum / 10)
  }
  if (carry != 0) carry.toString + sum else sum
}

def multByDigit(num: String, digit: Int) = {
  val (mult, carry) = num.foldRight(("", 0)) {
    case (d, (acc, carry)) =>
      val mult = d.asDigit * digit + carry
      ((mult % 10).toString + acc, mult / 10)
  }
  if (carry != 0) carry.toString + mult else mult
}

def mult(x: String, y: String) =
  y.foldLeft("")((acc, digit) => addNums(acc + "0", multByDigit(x, digit.asDigit)))

Sample:

scala> mult("18446744073709551616", "18446744073709551616")
res25: java.lang.String = 340282366920938463463374607431768211456
Works with: Scala version 2.8

Scala 2.8 introduces `scanLeft` and `scanRight` which can be used to simplify this further:

def adjustResult(result: IndexedSeq[Int]) = (
  result
  .map(_ % 10)        // remove carry from each digit
  .tail               // drop the seed carry
  .reverse            // put most significant digits on the left
  .dropWhile(_ == 0)  // remove leading zeroes
  .mkString
)

def addNums(x: String, y: String) = {
  val padSize = (x.length max y.length) + 1 // We want to keep a zero to the left, to catch the carry
  val paddedX = "0" * (padSize - x.length) + x
  val paddedY = "0" * (padSize - y.length) + y
  adjustResult((paddedX zip paddedY).scanRight(0) { 
    case ((dx, dy), last) => dx.asDigit + dy.asDigit + last / 10 
  })
}

def multByDigit(num: String, digit: Int) = adjustResult(("0"+num).scanRight(0)(_.asDigit * digit + _ / 10))

def mult(x: String, y: String) =
  y.foldLeft("")((acc, digit) => addNums(acc + "0", multByDigit(x, digit.asDigit)))

Scheme

Since Scheme already supports arbitrary precision arithmetic, build it out of church numerals. Don't try converting these to native integers. You will die waiting for the answer.

(define one (lambda (f) (lambda (x) (f x))))
(define (add a b) (lambda (f) (lambda (x) ((a f) ((b f) x)))))
(define (mult a b) (lambda (f) (lambda (x) ((a (b f)) x))))
(define (expo a b) (lambda (f) (lambda (x) (((b a) f) x))))
(define two (add one one))
(define six (add two (add two two)))
(define sixty-four (expo two six))
(display (mult (expo two sixty-four) (expo two sixty-four)))
Output:

(as run on Chicken Scheme on tio)

#<procedure (? f)>

Seed7

Seed7 supports arbitrary-precision arithmetic. The library bigint.s7i defines the type bigInteger. A bigInteger is a signed integer number of unlimited size. With library support the task can be solved by using the multiplication operator *:

$ include "seed7_05.s7i";
  include "bigint.s7i";

const proc: main is func
  begin
    writeln(2_**64 * 2_**64);
  end func;

Output:

340282366920938463463374607431768211456

This task seems to prefer an inferior implementation of a long multiplication, where long numbers are stored in decimal strings. Besides type safety there are seveal other drawbacks triggered by such a representation. E.g.: In almost all cases a representation with decimal strings leads to significant lower computing speed. The multiplication example below uses the requested inferior implementation:

$ include "seed7_05.s7i";

const func string: (in string: a) * (in string: b) is func
  result
    var string: product is "";
  local
    var integer: i is 1;
    var integer: j is 1;
    var integer: k is 0;
    var integer: carry is 0;
  begin
    if startsWith(a, "-") then
      if startsWith(b, "-") then
        product := a[2 ..] * b[2 ..];
      else
        product := "-" & a[2 ..] * b;
      end if;
    elsif startsWith(b, "-") then
      product := "-" & a * b[2 ..];
    else 
      product := "0" mult length(a) + length(b);
      for i range length(a) downto 1 do
        k := i + length(b);
        carry := 0; 
        for j range length(b) downto 1 do
          carry +:= (ord(a[i]) - ord('0')) * (ord(b[j]) - ord('0')) + (ord(product[k]) - ord('0'));
          product @:= [k] chr(carry rem 10 + ord('0'));
          carry := carry div 10;
          decr(k);
        end for;
        product @:= [k] chr(ord(product[k]) + carry);
      end for;
      while startsWith(product, "0") and length(product) >= 2 do
        product := product[2 ..];
      end while;
    end if;
  end func;

const proc: main is func
  begin
    writeln("-18446744073709551616" * "-18446744073709551616");
  end func;

The output is the same as with the superior solution.

Sidef

(Note that arbitrary precision arithmetic is native in Sidef).

say (2**64 * 2**64);
Translation of: Python
func add_with_carry(result, addend, addendpos) {
    loop {
        while (result.len < addendpos+1) {
            result.append(0)
        }
        var addend_digits = (addend.to_i + result[addendpos] -> to_s.chars)
        result[addendpos] = addend_digits.pop
        addend_digits.len > 0 || break
        addend = addend_digits.pop
        addendpos++
    }
}
 
func longhand_multiplication(multiplicand, multiplier) {
 
    var result = []
    var multiplicand_offset = 0
 
    multiplicand.reverse.each { |multiplicand_digit|
        var multiplier_offset = multiplicand_offset
        multiplier.reverse.each { |multiplier_digit|
            var multiplication_result = (multiplicand_digit.to_i * multiplier_digit.to_i -> to_s)
 
            var addend_offset = multiplier_offset
            multiplication_result.reverse.each { |result_digit_addend|
                add_with_carry(result, result_digit_addend, addend_offset)
                addend_offset++
            }
            multiplier_offset++
        }
        multiplicand_offset++
    }
 
    return result.join.reverse
}
 
say longhand_multiplication('18446744073709551616', '18446744073709551616')
Output:
340282366920938463463374607431768211456

Slate

(2 raisedTo: 64) * (2 raisedTo: 64).

Smalltalk

Note that arbitrary precision arithmetic is native in Smalltalk, and no-one would reinvent the wheel.

(2 raisedTo: 64) * (2 raisedTo: 64).

or, to display it:

Transcript showCR:(2 raisedTo: 64) * (2 raisedTo: 64).
"if ** is defined as alias: " Transcript showCR:(2 ** 64) * (2 ** 64).
Output:

340282366920938463463374607431768211456

There has been some discussion, if the above is fair... but then, I guess even using 32bit arithmetic is unfair to 8bit assembly language machines ;-) or comparing languages which have a 128bit integer against others... And is it fair, to change the challenge afterwards?

Anyway, here is a version which works on 2-digit decimal machine of the 1940's... (not that I know of any Smalltalk ever ported to a Zuse 1 :-)

"/ mhmh hard to avoid largeInteger arithmetic, 
"/ as the language does not specify, how many bits are used to represent
"/ SmallIntegers, and when the VM uses LargeInts.
"/ Lets assume, we run on a 2-digit decimal machine (smile).
"/ So lets work hard to avoid any convenient VM support,
"/ by doing decimal arithmetic (running on a decimal machine from the 1940s)
"/ and only allow 0..99 in a word (assuming it has a 2*2->4 digit multiply available)
"/ (smile: remember the Knuth MIX machine?)
"/ Long integers are represented as an array of such 2-digit words (least significant first).

"/ the code below should never ever been taken serious
"/ Not even as didactic example.
"/ NOONE WOULD EVER DO SUCH A STUPID THING

WORDMAX := 100.
add := 
    [:a :b |
        Array streamContents:[:s |
            |cy|
            cy := 0.
            1 to:(a size max:b size) do:[:wordIndex |
                |sum|
                wA := a at:wordIndex ifAbsent:0.
                wB := b at:wordIndex ifAbsent:0.
                sum := (wA + wB + cy).
                cy := (sum // WORDMAX).
                s nextPut:(sum % WORDMAX).
            ].
            cy ~~ 0 ifTrue:[s nextPut:cy].
        ].
    ].

"/ test 12,34 + 1
a := #( 34 12 ).
b := #( 1 ).
self assert:( add value:a value:b ) = #( 35 12 ).  

"/ test 99,99 + 1
a := #( 99 99 ).
b := #( 1  ).
self assert:( add value:a value:b ) = #( 00 00 1 ).  

"/ test 99,99,99,99 + 99,99,99,99
a := #( 99 99 99 99 ).
b := #( 99 99 99 99 ).
self assert:( add value:a value:b ) = #( 98 99 99 99 1).  

mulW := 
    [:a :w |
        |cy|
        cy := 0.
        Array streamContents:[:s |
            a do:[:wordA |
                |product|
                product := (wordA * w) + cy.
                s nextPut:(product % WORDMAX).
                cy := (product // WORDMAX)
            ].
            cy ~~ 0 ifTrue:[s nextPut:cy].
        ]
    ].

"/ test 1 * 2
a := #( 1 ).
self assert:( mulW value:a value:2) = #( 2).  

"/ test 2 * 99
a := #( 2 ).
self assert:( mulW value:a value:99) = #( 98 1).  

"/ test 99,99,99,99 * 99  
a := #( 99 99 99 99 ).
self assert:( mulW value:a value:99) = #( 01 99 99 99 98 ).  

mul := 
    [:a :b |
        |sum|

        sum := #( 0 ).
        b doWithIndex:[:wordB :wordIndex |
            partSum := mulW value:a value:wordB.
            shifted := (Array new:wordIndex-1 withAll:0),partSum.
            sum := add value:sum value:shifted.
        ].
        sum.
    ].

"/ test 99,99,99,99 * 99
a := #( 99 99 99 99 ).
b := #( 99  ).
self assert:( mul value:a value:b) = #( 01 99 99 99 98 ).

raise := 
    [:a :exp |
        |e rslt|

        rslt := #(1).
        t := a.
        e := exp.
        [e ~~ 0] whileTrue:[
            [(e bitAnd:1) == 0] whileTrue:[
                e := e bitShift:-1.
                t := mul value:t value:t.
            ].
            e := e - 1.
            rslt := mul value:rslt value:t.
        ].
        rslt.
    ].

"/ test 2 ** 64
a := #( 2  ).                                    
self assert:( raise value:a value:64) = #( 16 16 55 09 37 07 44 67 44 18).

"/ test (2 ** 64) * (2 ** 64)   
a := #( 2  ).                       
t := raise value:a value:64.
rslt := mul value:t value:t.           
self assert:rslt = #( 56 14 21 68 17 43 07 46 37 63 34 46 38 09 92 66 23 28 40 3).

"/ the biggest plus of having a decimal machine is that it makes printing so easy...
printOn := 
    [:n :stream |
        |first|
        first := true.
        n reverseDo:[:_2Digits |
            first 
                ifTrue:[ stream nextPutAll:(_2Digits printString)]    
                ifFalse:[ stream nextPutAll:(_2Digits printString leftPaddedTo:2 with:$0)]. 
            first := false.
        ].
    ].

printOn value:rslt value:Transcript.

"/ verify...
printedString := String streamContents:[:s | printOn value:rslt value:s].
self assert:(printedString = (2**64) squared printString)
Output:

3402823669293846346337467431768211456

The above code does not really integrate into the Smalltalk class library. For example, it will not allow mixed mode arithmetic between regular integers and Rosetta integers. Here is a full example in portable chunk file format which makes mixed mode arithmetic completely transparent (I implemented only addition and multiplication):

Integer
    subclass: #RosettaInteger
    instanceVariableNames:'digitArray'
    classVariableNames:'WORDMAX'
    package:'Rosetta demos'
!

!RosettaInteger class methodsFor:'initialization'!

initialize
    WORDMAX := 100.
! !

!RosettaInteger class methodsFor:'instance creation'!

newWithDigits:digitArray
    "returns a new RosettaInteger with a digitArray"

    ^ self basicNew digits:digitArray
!

fromInteger:anInteger
    "returns a new RosettaInteger with anInteger's value"

    |digits gen|

    gen := [:n :s |
            s nextPut:(n % 100).
            n > 99 ifTrue:[ gen value:(n // 100) value:s]].
    digits := Array streamContents:[:s | gen value:anInteger value:s].
    ^ self newWithDigits:digits
! !

!RosettaInteger class methodsFor:'helpers'!

addDigits:a and:b
    |add|

    add :=
        [:a :b |
            Array streamContents:[:s |
                |cy|
                cy := 0.
                1 to:(a size max:b size) do:[:wordIndex |
                    |sum|
                    wA := a at:wordIndex ifAbsent:0.
                    wB := b at:wordIndex ifAbsent:0.
                    sum := (wA + wB + cy).
                    cy := (sum // WORDMAX).
                    s nextPut:(sum % WORDMAX).
                ].
                cy ~~ 0 ifTrue:[s nextPut:cy].
            ].
        ].
    ^ add value:a value:b
!

mulDigits:a and:b
    |mulW|

    mulW :=
        [:a :w |
            |cy|
            cy := 0.
            Array streamContents:[:s |
                a do:[:wordA |
                    |product|
                    product := (wordA * w) + cy.
                    s nextPut:(product % WORDMAX).
                    cy := (product // WORDMAX)
                ].
                cy ~~ 0 ifTrue:[s nextPut:cy].
            ]
        ].

    mul :=
        [:a :b |
            |sum|

            sum := #( 0 ).
            b doWithIndex:[:wordB :wordIndex |
                partSum := mulW value:a value:wordB.
                shifted := (Array new:wordIndex-1 withAll:0),partSum.
                sum := self addDigits:sum and:shifted.
            ].
            sum.
        ].

    ^ mul value:a value:b.
! !

!RosettaInteger methodsFor:'private accessing'!
digits
    "return my digitArray"

    ^ digitArray
!

digits:digits
    "set my digitArray"

    digitArray := digits
! !

!RosettaInteger methodsFor:'arithmetic'!
+ aNumber
    ^ aNumber sumFromRosettaInteger:self
!

* aNumber
    ^ aNumber productFromRosettaInteger:self
!

raisedTo:exp
    |raise|

    raise :=
        [:a :exp |
            |e rslt|

            rslt := #(1).
            t := a.
            e := exp.
            [e ~~ 0] whileTrue:[
                [(e bitAnd:1) == 0] whileTrue:[
                    e := e bitShift:-1.
                    t := self class mulDigits:t and:t.
                ].
                e := e - 1.
                rslt := self class mulDigits:rslt and:t.
            ].
            rslt.
        ].
   ^ self class newWithDigits:(raise value:(self digits) value:exp)
!

sumFromRosettaInteger:anRInt
   ^ self class
        newWithDigits:(self class addDigits:(anRInt digits) and:(self digits))
!

productFromRosettaInteger:anRInt
   ^ self class newWithDigits:(self class mulDigits:(anRInt digits) and:(self digits))
! !

!RosettaInteger methodsFor:'printing'!

printOn:aStream
    |print|

    print :=
        [:n :stream |
            |first|
            first := true.
            n reverseDo:[:_2Digits |
                first
                    ifTrue:[ stream nextPutAll:(_2Digits printString)]
                    ifFalse:[ stream nextPutAll:(_2Digits printString leftPaddedTo:2 with:$0)].
                first := false.
            ].
        ].

    print value:(self digits) value:aStream
! !

!Integer methodsFor:'converting'!

asRosettaInteger
    ^ RosettaInteger fromInteger:self
! !

!Integer methodsFor:'double dispatching'!

sumFromRosettaInteger:anRInt
    ^ anRInt + (RosettaInteger fromInteger:self)
!

productFromRosettaInteger:anRInt
    ^ anRInt * (RosettaInteger fromInteger:self)
! !

RosettaInteger initialize
!

a := 124 asRosettaInteger.
e'a is: {a} ({a class})' printCR.
b := 333 asRosettaInteger.
e'b is: {b} ({b class})'printCR.
a_plus_b := a+b.
e'(a+b) is: {a_plus_b} ({(a_plus_b) class})' printCR.

c := 999 asRosettaInteger.
e'c is: {c} ({c class})' printCR.
c_plus_1 := c+1.
e'c+1 is: {c_plus_1} ({(c_plus_1) class})' printCR.

d := 100 asRosettaInteger.
e'd is: {d} ({d class})' printCR.
d_squared := d squared.
e'd squared is: {d_squared} ({d_squared class})' printCR.

e := 2 asRosettaInteger.
e_raisedTo_64 := e raisedTo:64.
e'2 raisedTo:64 is: {e_raisedTo_64} ({e_raisedTo_64 class})' printCR.

e_raisedTo_64_squared := (e raisedTo:64) squared.
e'result is: {e_raisedTo_64_squared} ({e_raisedTo_64_squared class})' printCR.

Transcript show:'once again: '. 
result := (2 asRosettaInteger raisedTo:64) squared.
Transcript showCR:result.
Output:
a is: 124 (RosettaInteger)
b is: 333 (RosettaInteger)
(a+b) is: 457 (RosettaInteger)
c is: 999 (RosettaInteger)
c+1 is: 1000 (RosettaInteger)
d is: 100 (RosettaInteger)
d squared is: 10000 (RosettaInteger)
2 raisedTo:64 is: 18446744073709551616 (RosettaInteger)
result is: 340282366920938463463374607431768211456 (RosettaInteger)
once again: 340282366920938463463374607431768211456

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:

package require Tcl 8.5

proc longmult {x y} {
    set digits [lreverse [split $x ""]]
    set result {0}
    set j -2
    foreach m [lreverse [split $y ""]] {
	set c 0
	set i [incr j]
	foreach d $digits {
	    set v [lindex $result [incr i]]
	    if {$v eq ""} {
		lappend result 0
		set v 0
	    }
	    regexp (.)(.)$ 0[expr {$v + $c + $d*$m}] -> c v
	    lset result $i $v
	}
	lappend result $c
    }
    # Reconvert digit list into a decimal number
    set result [string trimleft [join [lreverse $result] ""] 0]
    if {$result == ""} then {return 0} else {return $result}
}

puts [set n [expr {2**64}]]
puts [longmult $n $n]
puts [expr {$n * $n}]

outputs

18446744073709551616
340282366920938463463374607431768211456
340282366920938463463374607431768211456

uBasic/4tH

Translation of: Liberty BASIC
' now, count 2^64
Print "2^64"
a := "1"
For i = 1 To 64
  a = FUNC(_multByD (a, 2))
Next

Print Show (a)
Print "(check with known value)"
Print "18446744073709551616"
Print "(looks OK)"

' now let's do b*a stuff
Print
Print "2^64*2^64"
Print Show (FUNC (_longMult (a, a)))
Print "(check with known value)"
Print "340282366920938463463374607431768211456"
Print "(looks OK)"

End

' ---------------------------------------
_longMult
  Param (2)
  Local (7)

    c@ = 1
    If Peek (a@, 0) = Ord ("-") Then a@ = Chop (a@, 1) : c@ = -1
    d@ = 1
    If Peek (b@, 0) = Ord ("-") Then b@ = Chop (b@, 1) : d@ = -1

    e@ := ""
    f@ := ""
    g@ := ""
    For h@ = Len(a@)-1 To 0 Step -1
        i@ = Peek (a@, h@) - Ord ("0")
        f@ = FUNC(_multByD (b@, i@))
        e@ = FUNC(_addLong (e@, Join (f@, g@)))
        g@ = Join (g@, "0")
    Next                               ' print i@, f@, e@ 
                                       ' print e@    
    If c@*d@ < 0 Then e@ = Join ("-", e@)
Return (e@)

_multByD
  Param (2)                            ' multiply a@ by digit b@
  Local (5)
  
  c@ := ""
  d@ = 0

  For e@ = Len (a@)-1 To 0 Step -1
      f@ = Peek (a@, e@) - Ord ("0")
      g@ = f@*b@ + d@
      d@ = g@ / 10)
      g@ = g@ % 10
                                       ' print f@, g@
      c@ = Join (Str (g@), c@) 
  Next

  If d@ > 0 Then c@ = Join (Str (d@), c@)
Return (c@)                            ' print c@

_addLong                               ' add a@ + b@, for now only positive
  Param (2)
  Local (7)

    c@ = Max (Len(a@), Len(b@))
    a@ = FUNC(_Pad (a@, c@))
    b@ = FUNC(_Pad (b@, c@))
    
    d@ := ""                           ' result
    e@ = 0
    For f@ = c@-1 To 0 Step -1
        g@ = Peek (a@, f@) - Ord ("0")
        h@ = Peek (b@, f@) - Ord ("0")
        i@ = g@ + h@ + e@
        e@ = i@ / 10
        i@ = i@ % 10
                                       ' print g@, h@, i@
        d@ = Join (Str (i@), d@)
    Next
                                       ' print d@    
    If e@ > 0 Then d@ = Join (Str (e@), d@)
Return (d@)

_Pad                                   ' pad from right with 0 to length n
  Param (2)

  Do While Len (a@) < b@
    a@ = Join("0", a@)
  Loop
Return (a@)
Output:
2^64
18446744073709551616
(check with known value)
18446744073709551616
(looks OK)

2^64*2^64
340282366920938463463374607431768211456
(check with known value)
340282366920938463463374607431768211456
(looks OK)

0 OK, 0:478

UNIX Shell

In real shell scripts, I would use either `bc` or `dc` for this:

multiply() { echo "$1 $2 * p" | dc; }

But you can also do it with bash's built-in arithmetic:

add() { # arbitrary-precision addition
  local a="$1" b="$2" sum= carry=0
  if (( ${#a} < ${#b} )); then
    local t="$a" 
    a="$b" b="$t"
  fi

  while (( ${#a} )); do
    local -i d1="${a##${a%?}}" d2="10#0${b##${b%?}}" s=carry+d1+d2
    sum="${s##${s%?}}$sum"
    carry="10#0${s%?}"
    a="${a%?}" b="${b%?}"
  done
  echo "$sum"
}

multiply() { # arbitrary-precision multiplication
  local a="$1" b="$2" product=0
  if (( ${#a} < ${#b} )); then
    local t="$a" 
    a="$b" b="$t"
  fi

  local zeroes=
  while (( ${#b} )); do
    local m1="$a"
    local m2="${b##${b%?}}"
    local partial=$zeroes 
    local -i carry=0
    while (( ${#m1} )); do 
      local -i d="${m1##${m1%?}}"
      m1="${m1%?}"
      local -i p=d*m2+carry
      partial="${p##${p%?}}$partial"
      carry="10#0${p%?}"
    done
    partial="${carry#0}$partial"
    product="$(add "$product" "$partial")"
    zeroes=0$zeroes
    b="${b%?}"
  done
  echo "$product"
}
Output is the same either way:
$ multiply 18446744073709551616 18446744073709551616
340282366920928463463374607431768211456

Ursala

Natural numbers of unlimited size are a built in type, and arithmetic operations on them are available as library functions. However, since the task calls for explicitly implementing long multiplication, here is an implementation using nothing but language primitives. The numbers are represented as lists of booleans, LSB first. The compiler already knows how to parse and display them in decimal.

successor = ~&a^?\1! ~&ah?/~&NfatPRC ~&NNXatPC

sum = ~&B^?a\~&Y@a ~&B?abh/successor@alh2fabt2RC ~&Yabh2Ofabt2RC

product = ~&alrB^& sum@NfalrtPXPRCarh2alPNQX

x = 18446744073709551616

#show+

y = %nP product@iiX x

output:

340282366920938463463374607431768211456

Vedit macro language

This example multiplies the value on current line with the value on next line and stores result on the 3rd line.

BOL
#11 = EOL_Pos-Cur_Pos
#12 = EOL_Pos-1
Line(1)
#21 = EOL_Pos-Cur_Pos
#22 = EOL_Pos-1
EOL Ins_Newline
Ins_Char('0', COUNT, #11+#21)
#32 = Cur_Pos-1

for (#2 = 0; #2 < #21; #2++) {
    Goto_Pos(#22-#2) #5 = Cur_Char - '0'
    for (#1 = 0; #1 < #11; #1++) {
        Goto_Pos(#12-#1) #6 = Cur_Char - '0'
	#7 = #5 * #6
	#3 = #1 + #2
	while (#7 > 0) {
	    Goto_Pos(#32-#3)
	    #7 += Cur_Char - '0'
	    Ins_Char(#7%10 + '0', OVERWRITE)
	    #3++
	    #7 = #7/10 
        }
    }
}

Sample input and output:

18446744073709551616
18446744073709551616
0340282366920938463463374607431768211456

Visual Basic .NET

Translation of: C#

This uses the decimal type, (which has a MaxValue of 79,228,162,514,264,337,593,543,950,335). By limiting it to 10^28, it allows 28 decimal digits for the hi part, and 28 decimal digits for the lo part, 56 decimal digits total. A side computation of BigInteger assures that the results are accurate.

Imports System
Imports System.Console
Imports BI = System.Numerics.BigInteger

Module Module1

    Dim a As Decimal, mx As Decimal = 1E28D, hm As Decimal = 1E14D

    ' allows for 56 digit representation, using 28 decimal digits from each decimal
    Structure bd
        Public hi, lo As Decimal
    End Structure

    ' outputs bd structure as string, optionally inserting commas
    Function toStr(ByVal a As bd, ByVal Optional comma As Boolean = False) As String
        Dim r As String = If(a.hi = 0, String.Format("{0:0}", a.lo),
                                       String.Format("{0:0}{1:" & New String("0"c, 28) & "}", a.hi, a.lo))
        If Not comma Then Return r 
        Dim rc As String = ""
        For i As Integer = r.Length - 3 To 0 Step -3
            rc = "," & r.Substring(i, 3) & rc : Next
        toStr = r.Substring(0, r.Length Mod 3) & rc
            toStr = toStr.Substring(If(toStr.Chars(0) = "," , 1, 0))
    End Function

    ' needed because Math.Pow() returns a double
    Function Pow_dec(ByVal bas As Decimal, ByVal exp As UInteger) As Decimal
        If exp = 0 Then Pow_dec = 1D else Pow_dec = Pow_dec(bas, exp >> 1) : _
        Pow_dec *= Pow_dec : If (exp And 1) <> 0 Then Pow_dec *= bas
    End Function

    Sub Main(ByVal args As String())
         For p As UInteger = 64 To 95 - 1 Step 30                 ' show prescribed output and maximum power of 2 output
            Dim y As bd, x As bd : a = Pow_dec(2D, p)             ' init the bd variables, a = decimal value to be squared
            WriteLine("The square of (2^{0}):                    {1,38:n0}", p, a)
            x.hi = Math.Floor(a / hm) : x.lo = a Mod hm           ' setup for the squaring process
            Dim BS As BI = BI.Pow(CType(a, BI), 2)                ' for the BigInteger checking of result
            y.lo = x.lo * x.lo : y.hi = x.hi * x.hi               ' square the lo and the hi parts
            a = x.hi * x.lo * 2D                                  ' calculate the middle term (mid-term)
            y.hi += Math.Floor(a / hm) : y.lo += (a Mod hm) * hm  ' increment hi and lo parts with high and low parts of the mid-term
            While y.lo > mx : y.lo -= mx : y.hi += 1 : End While  ' check for overflow, adjust both parts as needed
            WriteLine(" is {0,75} (which {1} match the BigInteger computation)" & vbLf,
                toStr(y, True), If(BS.ToString() = toStr(y), "does", "fails to"))
        Next
    End Sub

End Module
Output:
Shown are the prescribed output and the maximum power of two that can be squared by this bd structure without overflowing.
The square of (2^64):                                18,446,744,073,709,551,616
 is                         340,282,366,920,938,463,463,374,607,431,768,211,456 (which does match the BigInteger computation)

The square of (2^94):                    19,807,040,628,566,084,398,385,987,584
 is 392,318,858,461,667,547,739,736,838,950,479,151,006,397,215,279,002,157,056 (which does match the BigInteger computation)

Wren

Translation of: Go
Library: Wren-fmt
import "./fmt" for Fmt

// argument validation
var d = Fn.new { |b|
    if (b < 48 || b > 57) Fiber.abort("digit 0-9 expected")
    return b - 48
}

// converts a list of bytes to a string
var b2s = Fn.new { |b| b.map { |e| String.fromByte(e) }.join() }

// add two numbers as strings
var add = Fn.new { |x, y|
    if (y.count > x.count) {
        var t = x
        x = y
        y = t
    }
    var b = List.filled(x.count+1, 0)
    var c = 0
    for (i in 1..x.count) {
        if (i <= y.count) c = c + d.call(y[y.count-i].bytes[0])
        var s = d.call(x[x.count-i].bytes[0]) + c
        c = (s/10).floor
        b[b.count-i] = (s%10) + 48
    }
    if (c == 0) return b2s.call(b[1..-1])
    b[0] = c + 48
    return b2s.call(b)
}

// multiply a number by a single digit
var mulDigit = Fn.new { |x, y|
    if (y == 48) return "0"
    y = d.call(y)
    var b = List.filled(x.count+1, 0)
    var c = 0
    for (i in 1..x.count) {
        var s = d.call(x[x.count-i].bytes[0]) * y + c
        c = (s/10).floor
        b[b.count-i] = (s%10) + 48
    }
    if (c == 0) return b2s.call(b[1..-1])
    b[0] = c + 48
    return b2s.call(b)
}

// multiply two numbers as strings
var mul = Fn.new { |x, y|
    var result = mulDigit.call(x, y[y.count-1].bytes[0])
    var zeros = ""
    var i = 2
    while (i <= y.count) {
        zeros = zeros + "0"
        result = add.call(result, mulDigit.call(x, y[y.count-i].bytes[0]) + zeros)
        i = i + 1
    }
    result = result.trimStart("0")
    if (result == "") result = "0"
    return result
}

var n = "18446744073709551616"
Fmt.print("$,s", mul.call(n, n))
Output:
340,282,366,920,938,463,463,374,607,431,768,211,456

XPL0

include c:\cxpl\stdlib;
char Two64, Product(40);
[Two64:= "18446744073709551616";
StrNMul(Two64, Two64, Product, 20);
Product(39):= Product(39)!$80;          \terminate string
Text(0, Product+1);                     \skip leading zero
]

Output:

340282366920938463463374607431768211456

zkl

[gnu] BigNums are supported via an extension library

var BN=Import("zklBigNum");
BN(2).pow(64) * BN(2).pow(64)
340282366920938463463374607431768211456

BN(2).pow(128) : "%,d".fmt(_)
340,282,366,920,938,463,463,374,607,431,768,211,456

   //42!, also BN(42).factorial()
[2..42].reduce(fcn(p,n){p*n},BN(1)) : "%,d".fmt(_)
1,405,006,117,752,879,898,543,142,606,244,511,569,936,384,000,000,000