|
| 1 | +/** |
| 2 | +* @license Apache-2.0 |
| 3 | +* |
| 4 | +* Copyright (c) 2026 The Stdlib Authors. |
| 5 | +* |
| 6 | +* Licensed under the Apache License, Version 2.0 (the "License"); |
| 7 | +* you may not use this file except in compliance with the License. |
| 8 | +* You may obtain a copy of the License at |
| 9 | +* |
| 10 | +* http://www.apache.org/licenses/LICENSE-2.0 |
| 11 | +* |
| 12 | +* Unless required by applicable law or agreed to in writing, software |
| 13 | +* distributed under the License is distributed on an "AS IS" BASIS, |
| 14 | +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 15 | +* See the License for the specific language governing permissions and |
| 16 | +* limitations under the License. |
| 17 | +*/ |
| 18 | + |
| 19 | +'use strict'; |
| 20 | + |
| 21 | +// VARIABLES // |
| 22 | + |
| 23 | +// Define a mask for the least significant 16 bits (low word): 65535 => 0x0000ffff => 00000000000000001111111111111111 |
| 24 | +var LOW_WORD_MASK = 0x0000ffff>>>0; // asm type annotation |
| 25 | + |
| 26 | + |
| 27 | +// MAIN // |
| 28 | + |
| 29 | +/** |
| 30 | +* Performs C-like multiplication of two unsigned 32-bit integers. |
| 31 | +* |
| 32 | +* ## Method |
| 33 | +* |
| 34 | +* - To emulate C-like multiplication without the aid of 64-bit integers, we recognize that a 32-bit integer can be split into two 16-bit words |
| 35 | +* |
| 36 | +* ```tex |
| 37 | +* a = w_h*2^{16} + w_l |
| 38 | +* ``` |
| 39 | +* |
| 40 | +* where \\( w_h \\) is the most significant 16 bits and \\( w_l \\) is the least significant 16 bits. For example, consider the maximum unsigned 32-bit integer \\( 2^{32}-1 \\) |
| 41 | +* |
| 42 | +* ```binarystring |
| 43 | +* 11111111111111111111111111111111 |
| 44 | +* ``` |
| 45 | +* |
| 46 | +* The 16-bit high word is then |
| 47 | +* |
| 48 | +* ```binarystring |
| 49 | +* 1111111111111111 |
| 50 | +* ``` |
| 51 | +* |
| 52 | +* and the 16-bit low word |
| 53 | +* |
| 54 | +* ```binarystring |
| 55 | +* 1111111111111111 |
| 56 | +* ``` |
| 57 | +* |
| 58 | +* If we cast the high word to 32-bit precision and multiply by \\( 2^{16} \\) (equivalent to a 16-bit left shift), then the bit sequence is |
| 59 | +* |
| 60 | +* ```binarystring |
| 61 | +* 11111111111111110000000000000000 |
| 62 | +* ``` |
| 63 | +* |
| 64 | +* Similarly, upon casting the low word to 32-bit precision, the bit sequence is |
| 65 | +* |
| 66 | +* ```binarystring |
| 67 | +* 00000000000000001111111111111111 |
| 68 | +* ``` |
| 69 | +* |
| 70 | +* From the rules of binary addition, we recognize that adding the two 32-bit values for the high and low words will return our original value \\( 2^{32}-1 \\). |
| 71 | +* |
| 72 | +* - Accordingly, the multiplication of two 32-bit integers can be expressed |
| 73 | +* |
| 74 | +* ```tex |
| 75 | +* \begin{align*} |
| 76 | +* a \cdot b &= ( a_h \cdot 2^{16} + a_l) \cdot ( b_h \cdot 2^{16} + b_l) \\ |
| 77 | +* &= a_l \cdot b_l + a_h \cdot b_l \cdot 2^{16} + a_l \cdot b_h \cdot 2^{16} + (a_h \cdot b_h) \cdot 2^{32} \\ |
| 78 | +* &= a_l \cdot b_l + (a_h \cdot b_l + a_l \cdot b_h) \cdot 2^{16} + (a_h \cdot b_h) \cdot 2^{32} |
| 79 | +* \end{align*} |
| 80 | +* ``` |
| 81 | +* |
| 82 | +* - We note that multiplying (dividing) an integer by \\( 2^n \\) is equivalent to performing a left (right) shift of \\( n \\) bits. |
| 83 | +* |
| 84 | +* - Further, as we want to return an integer of the same precision, for a 32-bit integer, the return value will be modulo \\( 2^{32} \\). Stated another way, we only care about the low word of a 64-bit result. |
| 85 | +* |
| 86 | +* - Accordingly, the last term, being evenly divisible by \\( 2^{32} \\), drops from the equation leaving the remaining two terms as the remainder. |
| 87 | +* |
| 88 | +* ```tex |
| 89 | +* a \cdot b = a_l \cdot b_l + (a_h \cdot b_l + a_l \cdot b_h) << 16 |
| 90 | +* ``` |
| 91 | +* |
| 92 | +* - Lastly, the second term in the above equation contributes to the middle bits and may cause the product to "overflow". However, we can disregard (`>>>0`) overflow bits due to modulo arithmetic, as discussed earlier with regard to the term involving the partial product of high words. |
| 93 | +* |
| 94 | +* @param {uinteger32} a - integer |
| 95 | +* @param {uinteger32} b - integer |
| 96 | +* @returns {uinteger32} product |
| 97 | +* |
| 98 | +* @example |
| 99 | +* var v = mul( 10>>>0, 4>>>0 ); |
| 100 | +* // returns 40 |
| 101 | +*/ |
| 102 | +function mul( a, b ) { |
| 103 | + var lbits; |
| 104 | + var mbits; |
| 105 | + var ha; |
| 106 | + var hb; |
| 107 | + var la; |
| 108 | + var lb; |
| 109 | + |
| 110 | + a >>>= 0; // asm type annotation |
| 111 | + b >>>= 0; // asm type annotation |
| 112 | + |
| 113 | + // Isolate the most significant 16-bits: |
| 114 | + ha = ( a>>>16 )>>>0; // asm type annotation |
| 115 | + hb = ( b>>>16 )>>>0; // asm type annotation |
| 116 | + |
| 117 | + // Isolate the least significant 16-bits: |
| 118 | + la = ( a&LOW_WORD_MASK )>>>0; // asm type annotation |
| 119 | + lb = ( b&LOW_WORD_MASK )>>>0; // asm type annotation |
| 120 | + |
| 121 | + // Compute partial sums: |
| 122 | + lbits = ( la*lb )>>>0; // asm type annotation; no integer overflow possible |
| 123 | + mbits = ( ((ha*lb) + (la*hb))<<16 )>>>0; // asm type annotation; possible integer overflow |
| 124 | + |
| 125 | + // The final `>>>0` converts the intermediate sum to an unsigned integer (possible integer overflow during sum): |
| 126 | + return ( lbits + mbits )>>>0; // asm type annotation |
| 127 | +} |
| 128 | + |
| 129 | + |
| 130 | +// EXPORTS // |
| 131 | + |
| 132 | +module.exports = mul; |
0 commit comments