Finite Field Multiplication

Question

Overview

Given the integer representation of three elements in GF(2^64), give the product of the first two elements over GF(2^64) with the reducing polynomial defined as the polynomial m such that m(x) = x^64 + n(x), where n is the polynomial representation of the third element.

You want me to what?

A finite field GF(p^q), where p is a prime and q is a natural number, contains p^q elements. (GF stands for Galois Field, which is another name for a finite field.) The elements of this finite field can be given many interpretations, but the two most common are as the integers between 0 and p^q-1, and as the polynomials with term coefficients between 0 and p-1 and degree less than q. There is a one-to-one correspondence between the base-p integer interpretations and the polynomial interpretations, from a one-to-one mapping between digits and term coefficients:

Integer interpretation (binary) | Polynomial interpretation (given x)
--------------------------------+------------------------------------
0                               | 0
1                               | 1
10                              | x
11                              | x + 1
100                             | x^2
101                             | x^2 + 1
1101                            | x^3 + x^2 + 1

Finite field arithmetic makes the most intuitive sense under the polynomial interpretation, but for brevity, the integer interpretation is used, and operations are described according to the integer representations of the finite field elements involved.

Addition in GF(p^q), like polynomial addition, operates with no carryover. Each digit in the sum is equal to the sum of digits in the addends, mod p. Likewise, subtraction in GF(p^q), like polynomial subtraction, is borrowless, with each digit in the difference being the difference in the corresponding digit in the minuend and subtrahend, mod p. If p = 2, both carryless addition and borrowless subtraction are equivalent to a bitwise XOR.

Multiplication in GF(p^q) relies on (but is not) polynomial multiplication, but with each term coefficient reduced mod p. If the factors are interpreted as integers rather than as polynomials, the operation can be implemented like long multiplication, but with summation of intermediate terms done without carry, as with finite field addition:

     1111
   * 1101
   ------
     1111
   1111
^ 1111
---------
  1001011

Multiplication in GF(p^q) also relies on a reducing polynomial, an irreducible polynomial with term coefficients in the set of integers between 0 and p-1, degree q, and leading coefficient 1 (and thus corresponding to an element in GF(p^(q+1))). Irreducible, in this context, means that it is not equal to any reduced polynomial product, as described above, of any other polynomials with term coefficients in the set of integers between 0 and p-1.

Last, multiplication in GF(p^q) relies on polynomial division, but with each term coefficient reduced mod p. Like reduced polynomial multiplication, this operation can be done directly on the integer representation by doing long division, but with borrowless subtraction rather than subtraction with borrow:

            1100
       ,--------
11111 / 10001001
      ^ 11111
      -------
         11100
       ^ 11111
       -------
            1101

Now, for the actual description of multiplication in GF(p^q).

First, the carryless product of the factors is calculated, as described above, keeping all digits in the product. Then, the remainder of the carryless division of this carryless product and the reducing polynomial (as interpreted as the integer interpretation of an element of GF(p^(q+1))) is taken. This remainder is the final product.

Task

Given 3 64-bit words A, B, and C, return pmod(pmul(A, B), C+1<<64)), where pmul is the carryless multiplication operation as described above, and pmod is the remainder of the carryless division operation as described above. C + 1<<64 can be assumed to correspond to an element in GF(2^65) that in turn represents an irreducible polynomial.

As with most code golf challenges, standard loopholes apply, and shortest program wins.

Test cases

Each test case is represented in four hexadecimal numbers, with the first three being the three arguments A, B, and C, as described above, and the fourth being the expected return value. (All values were calculated by Nemo, an open-source package for Julia that provides support for finite fields, among other things.)

0x0000000000000000 0x0000000000000000 0x000000000000001B 0x0000000000000000
0x0000000080000000 0x0000000080000000 0x000000000000001B 0x4000000000000000
0x0000000100000000 0x0000000100000000 0x000000000000001B 0x000000000000001B
0x0000000200000000 0x0000000200000000 0x000000000000001B 0x000000000000006C
0x0000000000000001 0x0000000000000001 0x702414004010004F 0x0000000000000002
0x0000000080000000 0x0000000100000000 0x702414004010004F 0x8000000000000000
0x0000000F00000000 0x0000000F00000000 0x702414004010004F 0x686D881284A01126
0x0000208000000000 0x0080000000000002 0x702414004010004F 0x870AF102748F8A93

Leaderboards

Here is a Stack Snippet to generate both a regular leaderboard and an overview of winners by language, copied from here:

var QUESTION_ID=160069,OVERRIDE_USER=8478;function answersUrl(e){return"https://api.stackexchange.com/2.2/questions/"+QUESTION_ID+"/answers?page="+e+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+ANSWER_FILTER}function commentUrl(e,s){return"https://api.stackexchange.com/2.2/answers/"+s.join(";")+"/comments?page="+e+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+COMMENT_FILTER}function getAnswers(){jQuery.ajax({url:answersUrl(answer_page++),method:"get",dataType:"jsonp",crossDomain:!0,success:function(e){answers.push.apply(answers,e.items),answers_hash=[],answer_ids=[],e.items.forEach(function(e){e.comments=[];var s=+e.share_link.match(/\d+/);answer_ids.push(s),answers_hash[s]=e}),e.has_more||(more_answers=!1),comment_page=1,getComments()}})}function getComments(){jQuery.ajax({url:commentUrl(comment_page++,answer_ids),method:"get",dataType:"jsonp",crossDomain:!0,success:function(e){e.items.forEach(function(e){e.owner.user_id===OVERRIDE_USER&&answers_hash[e.post_id].comments.push(e)}),e.has_more?getComments():more_answers?getAnswers():process()}})}function getAuthorName(e){return e.owner.display_name}function process(){var e=[];answers.forEach(function(s){var r=s.body;s.comments.forEach(function(e){OVERRIDE_REG.test(e.body)&&(r="<h1>"+e.body.replace(OVERRIDE_REG,"")+"</h1>")});var a=r.match(SCORE_REG);a&&e.push({user:getAuthorName(s),size:+a[2],language:a[1],link:s.share_link})}),e.sort(function(e,s){var r=e.size,a=s.size;return r-a});var s={},r=1,a=null,n=1;e.forEach(function(e){e.size!=a&&(n=r),a=e.size,++r;var t=jQuery("#answer-template").html();t=t.replace("{{PLACE}}",n+".").replace("{{NAME}}",e.user).replace("{{LANGUAGE}}",e.language).replace("{{SIZE}}",e.size).replace("{{LINK}}",e.link),t=jQuery(t),jQuery("#answers").append(t);var o=e.language;/<a/.test(o)&&(o=jQuery(o).text()),s[o]=s[o]||{lang:e.language,user:e.user,size:e.size,link:e.link}});var t=[];for(var o in s)s.hasOwnProperty(o)&&t.push(s[o]);t.sort(function(e,s){return e.lang>s.lang?1:e.lang<s.lang?-1:0});for(var c=0;c<t.length;++c){var i=jQuery("#language-template").html(),o=t[c];i=i.replace("{{LANGUAGE}}",o.lang).replace("{{NAME}}",o.user).replace("{{SIZE}}",o.size).replace("{{LINK}}",o.link),i=jQuery(i),jQuery("#languages").append(i)}}var ANSWER_FILTER="!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe",COMMENT_FILTER="!)Q2B_A2kjfAiU78X(md6BoYk",answers=[],answers_hash,answer_ids,answer_page=1,more_answers=!0,comment_page;getAnswers();var SCORE_REG=/<h\d>\s*([^\n,]*[^\s,]),.*?(\d+)(?=[^\n\d<>]*(?:<(?:s>[^\n<>]*<\/s>|[^\n<>]+>)[^\n\d<>]*)*<\/h\d>)/,OVERRIDE_REG=/^Override\s*header:\s*/i;

body{text-align:left!important}#answer-list,#language-list{padding:10px;width:290px;float:left}table thead{font-weight:700}table td{padding:5px}

<script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script> <link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b"> <div id="answer-list"> <h2>Leaderboard</h2> <table class="answer-list"> <thead> <tr><td></td><td>Author</td><td>Language</td><td>Size</td></tr></thead> <tbody id="answers"> </tbody> </table> </div><div id="language-list"> <h2>Winners by Language</h2> <table class="language-list"> <thead> <tr><td>Language</td><td>User</td><td>Score</td></tr></thead> <tbody id="languages"> </tbody> </table> </div><table style="display: none"> <tbody id="answer-template"> <tr><td>{{PLACE}}</td><td>{{NAME}}</td><td>{{LANGUAGE}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr></tbody> </table> <table style="display: none"> <tbody id="language-template"> <tr><td>{{LANGUAGE}}</td><td>{{NAME}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr></tbody> </table>

Answer 1

APL (Dyalog Unicode), ³⁹ 36 bytes^SBCS

{×≢⍵:(⍺×⊃⌽⍵)≠((⍺⍺×⊃⍺)≠1↓⍺,0)∇¯1↓⍵⋄0}

Try it online!

^{-3 bytes thanks to @ngn.}

A recursive inline operator (dop). You can assign it to f and use it as a(m f)b.

A port of orlp's Python answer, but using bit arrays instead of integers. Dyalog APL has neither 64-bit integers nor bitwise ops, but can simulate them pretty well using array ops.

How it works

{×≢⍵:(⍺×⊃⌽⍵)≠((⍺⍺×⊃⍺)≠1↓⍺,0)∇¯1↓⍵⋄0}  ⍝ ⍺←A, ⍵←B, ⍺⍺←M
{×≢⍵:                            ⋄ }  ⍝ If B has at least one bit,
                            ∇         ⍝   Recursive call with M fixed with...
                      1↓⍺,0           ⍝     A shifted left once
                     ≠                ⍝     XORed with
              (⍺⍺×⊃⍺)                 ⍝     M × highest bit of A
             (             )          ⍝     as left arg, and
                             ¯1↓⍵     ⍝     B shifted right once as right arg
     (⍺×⊃⌽⍵)≠                         ⍝   XOR the result with A × LSB of B
                                  0   ⍝ Otherwise, return a zero,
                                      ⍝ which is equivalent to 64 zero bits

The single zero is returned only on recursion. Then it becomes the right arg of (⍺×⊃⌽⍵)≠, and it acts exactly like a vector of 64 zeros (since (⍺×⊃⌽⍵) evaluates to a 64-length vector).

Answer 2

Perl 5 -p, 51 bytes

$\=$'*($\>>63)^$\<<1^$`*($&<<$_&1)for-63..!/ .* /}{

Try it online!

I originally thought this would be hard for perl, but this isn't too bad. All the $s of course still hurt for arithmetic challenges..

Answer 3

Python 2, 57 bytes

f=lambda a,b,m:a*b and b%2*a^f(a*2%2**64^(a>>63)*m,b/2,m)

I really tried to golf this in Rust, but the language is just not suitable for golfing. So here it is in Python.

Answer 4

Jelly, 22 bytes

Bæc/ḣ-64æcB}ʋƬUSḣ64UḂḄ

Try it online!

Works either as a function or a full program. This takes two arguments: a 2-element list [A,B] for the numbers to multiply; and the argument C that specifies the polynomial. The TIO link contains a footer that prints the result in hexadecimal, making it easier to verify the test case.

Algorithm

The calculation is done by representing each finite field element that we're working on as a list of bits that represent the coefficients of the polynomial. (For working on finite fields with this representation, we'd really want them to be little-endian, but Jelly's base conversion operators are big-endian, and thus the program ends up swapping endianness a couple of times.) Converting A and B to this form, then convolving them (which is the same operation as polynomial multiplication), gives us a polynomial (not necessarily a normalised one) representing the result.

We then need to normalize it (i.e. reduce it to degree at most 63) by doing the polynomial long division, i.e. taking the result modulo x⁶⁴+C (where C is interpreted as a polynomial). This is the same as taking the result modulo x⁶⁴−C because addition and subtraction are equivalent in binary finite fields. If we right-shift the result polynomial by 64 degrees, then multiply by C and add back the parts of the polynomial that fell off during the shift, this doesn't change the value modulo x⁶⁴−C (because some multiple of x⁶⁴ got subtracted from the result and the same multiple of C got added, thus reducing the result by a multiple of x⁶⁴−C and thus not changing the remainder upon dividing by x⁶⁴−C). Additionally, because C has degree at most 63 (it's a "64-bit word" interpreted as a polynomial), the shift-multiply-add sequence must reduce the degree of the resulting polynomial unless it's already 63 or lower, i.e. it's already normalized. It's therefore possible to repeat the shift-multiply-add (with polynomial addition and multiplication in GF(2), i.e. pointwise XOR and XOR-convolution) until we reach a value that remains unchanged, which must be the normalized result.

This algorithm is tweaked a little to be terser in Jelly. There are only two major changes. One is that we don't bother using only 0 and 1 as coefficients in the polynomials; in GF(2), 0 and 2 are the same number (all that matters is oddness versus evenness), so we just use regular pointwise addition (which doesn't carry from one coefficient into the next) instead of pointwise XOR and regular convolution instead of XOR-convolution (this saves lots of bytes during the program and requires only a Ḃ to fix the coefficients back to the 0-or-1 representation near the end). The other is that instead of immediately adding back the coefficients that got shifted off, they're remembered in a list and added back at the end (which is clearly equivalent because if we added them back immediately, they would just get shifted off again on the next iteration of the shift-multiply-add loop). (We actually remember the whole polynomial, but just take the least significant 64 coefficients for summing purposes.)

In order for the sum to work correctly, we need to use a little-endian representation (so that the coefficients that are added to each other when summing, i.e. the ones that are the same distance from the left end of the list, correspond to the same power of x0). It's possible to do the sum and the cut down to 64 coefficients either way round, and doing the sum first saves on the need for an extra quick to specify which array rank we want to work on. Then we convert the polynomial back to a big-endian representation, and then to a number, as the program's output.

Explanation

Bæc/ḣ-64æcB}ʋƬUSḣ64UḂḄ
B                        Convert {[A,B]} to big-endian binary
 æc                      Take the convolution of
   /                       the converted version of A with that of B
             Ƭ           Repeatedly loop, until nothing changes:
            ʋ              (grouping/bracketing mark, sets the Ƭ's scope)
    ḣ                      take the prefix
     -                       that omits
      64                     the last 64 elements
        æc                 and take its convolution with
          B                  the big-endian binary representation of
           }                   the second argument
              U          Reverse {the result of} each {iteration}
               S         Sum corresponding elements {of the reversed results}
                ḣ64      Take the first 64 elements of the sum
                   U     Reverse back
                    Ḃ    Replace each element with its value mod 2
                     Ḅ   Convert from big-endian binary to decimal

The Us mark where the program changes from big-endian to little-endian notation for its polynomial, and where it changes back again.

Answer 5

MMIX, 36 bytes (9 instrs)

00000000: e3030000 76040100 3f010101 70ff0003  ẉ¤¡¡v¥¢¡?¢¢¢p”¡¤
00000010: 3b000001 c6030304 c60000ff 5b01fffa  ;¡¡¢İ¤¤¥İ¡¡”[¢”«
00000020: f8040000                             ẏ¥¡¡

Disassembled and commented

ffm SET  $3,0           // r = 0
0H  ZSOD $4,$1,$0       // loop: t0 = B & 1? A : 0
    SRU  $1,$1,1        // B >>= 1
    ZSN  $255,$0,$3     // t1 = A >> 63? C : 0
    SLU  $0,$0,1        // A <<= 1
    XOR  $3,$3,$4       // r ^= t0
    XOR  $0,$0,$255     // A ^= t1
    PBNZ $1,0B          // if(B) goto loop
    POP  4,0            // return r

Russian peasant multiplication. Surprisingly simple! I think this is the first time MMIX has tied for shortest.

Answer 6

PARI/GP, 69 bytes

f(a,b,c)=lift(g(a)*g(b)%(x^64+g(c)))%(x-2)
g(n)=Pol(Mod(binary(n),2))

Attempt This Online!