x86-64 machine code, 34 bytes
Calling convention = x86-64 System V x32 ABI (register args with 32-bit pointers in long mode).
The function signature is void stewie_x87_1reg(float *seq_buf, unsigned Nterms);. The function receives the x0 and x1 seed values in the first two elements of the array, and extends the sequence out to at least N more elements. The buffer must be padded out to 2 + N-rounded-up-to-next-multiple-of-4. (i.e. 2 + ((N+3)&~3), or just N+5).
Requiring padded buffers is normal in assembly for high-performance or SIMD-vectorized functions, and this unrolled loop is similar, so I don't think it's bending the rules too far. The caller can easily (and should) ignore all the padding elements.
Passing x0 and x1 as a function arg not already in the buffer would cost us only 3 bytes (for a movlps [rdi], xmm0 or movups [rdi], xmm0), although this would be a non-standard calling convention since System V passes struct{ float x,y; }; in two separate XMM registers.
This is objdump -drw -Mintel output with a bit of formatting to add comments
0000000000000100 <stewie_x87_1reg>:
;; load inside the loop to match FSTP at the end of every iteration
;; x[i-1] is always in ST0
;; x[i-2] is re-loaded from memory
100: d9 47 04 fld DWORD PTR [rdi+0x4]
103: d8 07 fadd DWORD PTR [rdi]
105: d9 57 08 fst DWORD PTR [rdi+0x8]
108: 83 c7 10 add edi,0x10 ; 32-bit pointers save a REX prefix here
10b: d8 4f f4 fmul DWORD PTR [rdi-0xc]
10e: d9 57 fc fst DWORD PTR [rdi-0x4]
111: d8 6f f8 fsubr DWORD PTR [rdi-0x8]
114: d9 17 fst DWORD PTR [rdi]
116: d8 7f fc fdivr DWORD PTR [rdi-0x4]
119: d9 5f 04 fstp DWORD PTR [rdi+0x4]
11c: 83 ee 04 sub esi,0x4
11f: 7f df jg 100 <stewie_x87_1reg>
121: c3 ret
0000000000000122 <stewie_x87_1reg.end>:
## 0x22 = 34 bytes
This C reference implementation compiles (with gcc -Os) to somewhat similar code. gcc picks the same strategy I did, of keeping just one previous value in a register.
void stewie_ref(float *seq, unsigned Nterms)
{
for(unsigned i = 2 ; i<Nterms ; ) {
seq[i] = seq[i-2] + seq[i-1]; i++;
seq[i] = seq[i-2] * seq[i-1]; i++;
seq[i] = seq[i-2] - seq[i-1]; i++;
seq[i] = seq[i-2] / seq[i-1]; i++;
}
}
I did experiment with other ways, including a two-register x87 version that has code like:
; part of loop body from untested 2-register version. faster but slightly larger :/
; x87 FPU register stack ; x1, x2 (1-based notation)
fadd st0, st1 ; x87 = x3, x2
fst dword [rdi+8 - 16] ; x87 = x3, x2
fmul st1, st0 ; x87 = x3, x4
fld st1 ; x87 = x4, x3, x4
fstp dword [rdi+12 - 16] ; x87 = x3, x4
; and similar for the fsubr and fdivr, needing one fld st1
You'd do it this way if you were going for speed (and SSE wasn't available)
Putting the loads from memory inside the loop instead of once on entry might have helped, since we could just store the sub and div results out of order, but still it needs two FLD instructions to set up the stack on entry.
I also tried using SSE/AVX scalar math (starting with values in xmm0 and xmm1), but the larger instruction size is killer. Using addps (since that's 1B shorter than addss) helps a tiny bit. I used AVX VEX-prefixes for non-commutative instructions, since VSUBSS is only one byte longer than SUBPS (and the same length as SUBSS).
; untested. Bigger than x87 version, and can spuriously raise FP exceptions from garbage in high elements
addps xmm0, xmm1 ; x3
movups [rdi+8 - 16], xmm0
mulps xmm1, xmm0 ; xmm1 = x4, xmm0 = x3
movups [rdi+12 - 16], xmm1
vsubss xmm0, xmm1, xmm0 ; not commutative. Could use a value from memory
movups [rdi+16 - 16], xmm0
vdivss xmm1, xmm0, xmm1 ; not commutative
movups [rdi+20 - 16], xmm1
Tested with this test-harness:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
int main(int argc, char**argv)
{
unsigned seqlen = 100;
if (argc>1)
seqlen = atoi(argv[1]);
float first = 1.0f, second = 2.1f;
if (argc>2)
first = atof(argv[2]);
if (argc>3)
second = atof(argv[3]);
float *seqbuf = malloc(seqlen+8); // not on the stack, needs to be in the low32
seqbuf[0] = first;
seqbuf[1] = second;
for(unsigned i=seqlen ; i<seqlen+8; ++i)
seqbuf[i] = NAN;
stewie_x87_1reg(seqbuf, seqlen);
// stewie_ref(seqbuf, seqlen);
for (unsigned i=0 ; i< (2 + ((seqlen+3)&~3) + 4) ; i++) {
printf("%4d: %g\n", i, seqbuf[i]);
}
return 0;
}
Compile with nasm -felfx32 -Worphan-labels -gdwarf2 golf-stewie-sequence.asm &&
gcc -mx32 -o stewie -Og -g golf-stewie-sequence.c golf-stewie-sequence.o
Run the first test-case with ./stewie 8 1 3
If you don't have x32 libraries installed, use nasm -felf64 and leave gcc using the default -m64. I used malloc instead of float seqbuf[seqlen+8] (on the stack) to get a low address without having to actually build as x32.
Fun fact: YASM has a bug: it uses a rel32 jcc for the loop branch, when the branch target has the same address as a global symbol.
global stewie_x87_1reg
stewie_x87_1reg:
;; ended up moving all prologue code into the loop, so there's nothing here
.loop:
...
sub esi, 4
jg .loop
assembles to ... 11f: 0f 8f db ff ff ff jg 100 <stewie_x87_1reg>
Nbe 0-based? So take as inputs 1 less than the N shown in your examples. I guess taking N-2 is too much to ask for... :-P \$\endgroup\$