Terje Mathisen wrote:

>
> Something like

> if (n < 11) return ((0<<60)+(1<<54)+(1<<48)+(2<<42)...(21<<6)+34)>>
> (60-n*6) & 63;

unsigned char bytes[14] = { 0, 1, 1, 2, 3, 5, 8, 13, 34, 55, 89, 144, 233 };
unsigned short half[11] = {377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368 };

fib( unsigned n )
{ if( n <= 14 ) return bytes[n];
if( n <= 25 ) return half[n-14];
return fib( n-1 ) + fib( n-2 );
}

MitchAlsup wrote:
> Terje Mathisen wrote:
>
>>
>> Something like
>
>>    if (n < 11) return ((0<<60)+(1<<54)+(1<<48)+(2<<42)...(21<<6)+34)>>
>> (60-n*6) & 63;
>
>
> unsigned char bytes[14] = { 0, 1, 1, 2, 3, 5, 8, 13, 34, 55, 89, 144,
> 233 };
> unsigned short half[11] = {377, 610, 987, 1597, 2584, 4181, 6765, 10946,
> 17711, 28657, 46368 };
>
> fib( unsigned n )
> {
>     if( n <= 14 ) return bytes[n];
>     if( n <= 25 ) return half[n-14];
>     return fib( n-1 ) + fib( n-2 );
> }
Very nice except for not using the unroll-by-two, started with the two
last table entries? (Or second-to-last to fix any parity issue)

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

unsigned char bytes[14] = { 0, 1, 1, 2, 3, 5, 8, 13, 34, 55, 89, 144, 233 };
unsigned short half[11] = {377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368 };

fib( unsigned n )
{ if( n <= 14 ) return bytes[n];
if( n <= 25 ) return half[n-14];
uint64_t a = half[9], b = half[10];
for( n -= 25; n; n-- )
{
b += a;
a += b;
}
return a;
}

mitchalsup@aol.com (MitchAlsup) writes:

> Tim Rentsch wrote:
>
>> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>
>>>
>
>>>> Incidentally your code produces wrong values, as
>>>> fib(0) is 0 and fib(1) is 1.
>>>
>>> So return n for (n<2)?
>>
>> In situations like this one, usually I would favor something
>> more like this:
>>
>> if( n < 8 ) return 0xd8532110u >> n*4 & 15;
>
> Which is::
>
> CMP R2,R1,#8
> PLT R2,TTT
> LEA R2,[IP,R1<<2,0xF00000000-.]
> SLL R1,0xd853110,R2
> RET
>
> ....

Probably not exactly that, but thank you.

Terje Mathisen <terje.mathisen@tmsw.no> writes:

> Tim Rentsch wrote:
>
>> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>>
>>> Tim Rentsch wrote:
>>>
>>>> Thomas Koenig <tkoenig@netcologne.de> writes:
>>
>> [...]
>>
>>>>> unsigned long int fib (int n)
>>>>> {
>>>>> unsigned long f1, f2, fn;
>>>>>
>>>>> if (n < 2)
>>>>> return 1;
>>>>>
>>>>> f1 = 1;
>>>>> f2 = 1;
>>>>> while (n >= 2) {
>>>>> fn = f2 + f1;
>>>>> f1 = f2;
>>>>> f2 = fn;
>>>>> n--;
>>>>> }
>>>>> return fn;
>>>>> }
>>
>> [...]
>>
>>>> Incidentally your code produces wrong values, as
>>>> fib(0) is 0 and fib(1) is 1.
>>>
>>> So return n for (n<2)?
>>
>> In situations like this one, usually I would favor something
>> more like this:
>>
>> if( n < 8 ) return 0xd8532110u >> n*4 & 15;
>
> Nice!
>
> Could be extended a little bit more with a lookup table, or to 16
> values with a SIMD in-register lookup operation.

No reason to provide more values. It turns out these initial 8
values are just what the doctor ordered as a lead in to my latest
iterative version. In earlier versions putting in this special
case hurt more than it helped.

> Using a 64-bit
> constant allows 11 values to fit inside 6-bit sub-containers, as long
> as they are stored in reverse order so that the bottom value (0) can
> survive without the top two bits.
>
> Something like
>
> if (n < 11) return ((0<<60)+(1<<54)+(1<<48)+(2<<42)...(21<<6)+34)>>
> (60-n*6) & 63;

Probably more clever than useful. If I were doing this most likely
I would write the constant in octal:

0001010203051015254267u >> (60-6*n) & 63

Note by the way that you left off the value 55, which puts the
values 21 and 34 in the wrong places.

> Anyway, this is code golfing and not a realistic scenario where you
> would definitely use a direct lookup table for all the possible
> values.

A full lookup table is too doggone big. Also part of the point
of using the constant shift method is not to access data memory,
as that has other costs (and which are unlikely to show up in
speed measurements). The computational form is plenty fast
enough (worst case time between 7 and 8 nanoseconds in my tests).

This technique of packing several small values into a constant
and shifting to get a specific result is useful in lots of
situations. For example, if we are writing a function to test
for primality, it can start with

if( n < 64 ) return 0240404242024042424254 >>n &1;

for a fast test on small integers.

mitchalsup@aol.com (MitchAlsup) writes:

> Terje Mathisen wrote:
>
>> Something like
>>
>> if (n < 11) return
>> ((0<<60)+(1<<54)+(1<<48)+(2<<42)...(21<<6)+34)>> (60-n*6) & 63;
>
> unsigned char bytes[14] = { 0, 1, 1, 2, 3, 5, 8, 13, 34, 55, 89, 144, 233 };
> unsigned short half[11] = {377, 610, 987, 1597, 2584, 4181, 6765,
> 10946, 17711, 28657, 46368 };
>
> fib( unsigned n )
> {
> if( n <= 14 ) return bytes[n];
> if( n <= 25 ) return half[n-14];
> return fib( n-1 ) + fib( n-2 );
> }

Several bugs in this code (and exponential blowup when n > 25).

mitchalsup@aol.com (MitchAlsup) writes:

> unsigned char bytes[14] = { 0, 1, 1, 2, 3, 5, 8, 13, 34, 55, 89, 144, 233 };
> unsigned short half[11] = {377, 610, 987, 1597, 2584, 4181, 6765,
> 10946, 17711, 28657, 46368 };
>
> fib( unsigned n )
> {
> if( n <= 14 ) return bytes[n];
> if( n <= 25 ) return half[n-14];
> uint64_t a = half[9], b = half[10];
> for( n -= 25; n; n-- )
> {
> b += a;
> a += b;
> }
> return a;
> }

You need to divide n-25 by 2, and adjust for an even or odd
value before the divide.

Tim Rentsch wrote:
> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>
>> Tim Rentsch wrote:
>>
>>> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>>>
>>>> Tim Rentsch wrote:
>>>>
>>>>> Thomas Koenig <tkoenig@netcologne.de> writes:
>>>
>>> [...]
>>>
>>>>>> unsigned long int fib (int n)
>>>>>> {
>>>>>> unsigned long f1, f2, fn;
>>>>>>
>>>>>> if (n < 2)
>>>>>> return 1;
>>>>>>
>>>>>> f1 = 1;
>>>>>> f2 = 1;
>>>>>> while (n >= 2) {
>>>>>> fn = f2 + f1;
>>>>>> f1 = f2;
>>>>>> f2 = fn;
>>>>>> n--;
>>>>>> }
>>>>>> return fn;
>>>>>> }
>>>
>>> [...]
>>>
>>>>> Incidentally your code produces wrong values, as
>>>>> fib(0) is 0 and fib(1) is 1.
>>>>
>>>> So return n for (n<2)?
>>>
>>> In situations like this one, usually I would favor something
>>> more like this:
>>>
>>> if( n < 8 ) return 0xd8532110u >> n*4 & 15;
>>
>> Nice!
>>
>> Could be extended a little bit more with a lookup table, or to 16
>> values with a SIMD in-register lookup operation.
>
> No reason to provide more values. It turns out these initial 8
> values are just what the doctor ordered as a lead in to my latest
> iterative version. In earlier versions putting in this special
> case hurt more than it helped.
>
>> Using a 64-bit
>> constant allows 11 values to fit inside 6-bit sub-containers, as long
>> as they are stored in reverse order so that the bottom value (0) can
>> survive without the top two bits.
>>
>> Something like
>>
>> if (n < 11) return ((0<<60)+(1<<54)+(1<<48)+(2<<42)...(21<<6)+34)>>
>> (60-n*6) & 63;
>
> Probably more clever than useful. If I were doing this most likely
> I would write the constant in octal:
>
> 0001010203051015254267u >> (60-6*n) & 63
>
> Note by the way that you left off the value 55, which puts the
> values 21 and 34 in the wrong places.
>
>> Anyway, this is code golfing and not a realistic scenario where you
>> would definitely use a direct lookup table for all the possible
>> values.
>
> A full lookup table is too doggone big. Also part of the point
> of using the constant shift method is not to access data memory,
> as that has other costs (and which are unlikely to show up in
> speed measurements). The computational form is plenty fast
> enough (worst case time between 7 and 8 nanoseconds in my tests).

I agree, this is indeed "fast enough".
>
> This technique of packing several small values into a constant
> and shifting to get a specific result is useful in lots of
> situations. For example, if we are writing a function to test
> for primality, it can start with
>
> if( n < 64 ) return 0240404242024042424254 >>n &1;
>
> for a fast test on small integers.
>

I have been using this particular trick since at least the 386, so
around 1990. :-)

It did not make the same kind of sense when registers only had 16 bits
and a table lookup was about the same cost as a shift & mask.

Your 64-bit bitmap here is (inversely) related to the modulo 30 packing
I'm using in my sieve, where the 8 possible prime locations in each
block of 30 numbers will fit in a byte.

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

MitchAlsup wrote:
> unsigned char bytes[14] = { 0, 1, 1, 2, 3, 5, 8, 13, 34, 55, 89, 144,
> 233 };
> unsigned short half[11] = {377, 610, 987, 1597, 2584, 4181, 6765, 10946,
> 17711, 28657, 46368 };
>
> fib( unsigned n )
> {
>     if( n <= 14 ) return bytes[n];
>     if( n <= 25 ) return half[n-14];
>     uint64_t a = half[9], b = half[10];
>     for( n -= 25; n; n-- )
>     {
>          b += a;
>          a += b;
>     }
>     return a;
> }

As Tim wrote, you need something like

uint64_t a = half[8+(n&1)], b = half[9+(n&1)];
for (n = (n-24)>>1; n; n--) {}

to start the unrolled loop with the correct parity, although I might
have an off-by-one error in the (a,b) initializer here...

Since Tim is using a 4-way unroll, I believe he would be loading 4
values from a table of at least 7 entries to handle loop alignement.

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

Terje Mathisen wrote:
>
> Since Tim is using a 4-way unroll, I believe he would be loading 4
> values from a table of at least 7 entries to handle loop alignement.
>

I took another look at this, and it looks pretty nice to use 4 starting
entries (a..d in increasing order), then a two-cycle recurrence:

a = c+d
b = c+2*d
c = a+b
d = a+2*b

which every LEA-capable compiler around could convert into

lea rax,[rcx+rdx]
lea rbx,[rcx+2*rdx]

lea rcx,[rax+rbx]
lea rdx,[rax+2*rbx]

The pairs of instructions can run at the same time, so this can run in
two cycles as long as back-to-back LEA can run in subsequent cycles.

To run even faster than this I seem to need actual multiplications, with
small fibonacci numbers as the multipliers:

f(n+1) = f(n)+f(n-1)
f(n+2) = 2*f(n)+f(n-1)
f(n+3) = 3*f(n)+2*f(n-1)
f(n+4) = 5*f(n)+3*(f(n-1)
....
f(n+8) = 34*(f(n)+21*f(n-1)

Terje
--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

Tim Rentsch wrote:
> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>
>> Tim Rentsch wrote:
>>
>>> Daniel James <daniel@me.invalid> writes:
>>>
>>>> On 31/12/2023 09:12, Pancho wrote:
>>>>
>>>>> In fact, given that fib grows exponentially, so we only ever need to
>>>>> calculate it for small values of n, the only ?bad? solution is the
>>>>> naive double recursive one.
>>>>
>>>> Well, I might disagree.
>>>>
>>>> Given that we only ever need to calculate it for small values of n (as
>>>> you say, and if that's true) its efficiency doesn't matter. If you
>>>> need to compute it often enough that efficiency might become relevant
>>>> just put the results in a table and use that. Given that fib goes over
>>>> 10^12 after about 60 terms it needn't be a very big table.
>>>
>>> When the size of the table is more than ten times the size of
>>> the code, and the code is reasonably fast, using a table is a
>>> poor choice. Fibonacci functions don't have to be slow.
>>>
>>>> The naive doubly recursive solution has the huge merit of being
>>>> human-readable.
>>>
>>> But unacceptably slow.
>>>
>>>> Unlike, say:
>>>>
>>>> double fib( int n )
>>>> {
>>>> return (pow( (1+sqrt(5.0))/2, n+1 )
>>>> - pow( (1-sqrt(5.0))/2, n+1 ))
>>>> / sqrt(5.0);
>>>> }
>>>>
>>>> ... which is great for fibs of big values, but not fast for small ones.
>>>
>>> Unacceptable for fibs of big values, because the values computed
>>> are wrong.
>>
>> So either you move to quad or arbitrary precision (which will of
>> course be somewhat slow,
>
> To me, what makes floating point unattractive in this case is
> great care must be taken to ensure correct values are produced.
> Add to that the increased time cost, and using floating point
> just looks like a lot more downside than upside.
>
>> or you use a hybrid technique, with lookup
>> tables for smallish n, iterative loop for somewhat larger and your
>> O(log(n)) approach only for very large n?
>
> Oh, you remembered I wrote an O(log(n)) fibonacci function. That
> is nice. :)
>
> Based on what I know now, my inclination is to use the constant-shift
> technique (as seen in another post) for n < 8, the O(log(n)) method
> for results that will be more than 64 bits, and a small-constant-linear
> method in between those two regions. (Incidentally, the intermediate
> method is something I had already written, and might be seen as a
> generalization of your two-fold speedup code. You might want to take
> a look at enhancing the two-fold speedup idea so it does at most one of
> those, and then after that gets a four-fold speedup. The best
> performing code I have for 64 bits goes two steps farther, so all the
> iterations after the "startup transient" are 16-fold speedup.)

Do you actually get 16 results per clock cycle?

I don't see _any_ way to do that...
>
>> Working with 256-bit uints, each addition will take 4-8 clock cycles,
>> and it overflows well before n reaches 400, right? In that case a
>> worst case iterative solution is less than 3000 clock cycles or about
>> a microsecond away.
>
> By test, my O(log(n)) code runs in 25-30 nanoseconds for 128-bit

That is impressive!

> results. At a guess, 256-bit results would run about 5 to 10
> times slower than that.
>
> By the way, your two-fold speedup code runs about 1.6 times as
> fast as using simple iteration.
>
That's in the expected ballpark.

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

Terje Mathisen <terje.mathisen@tmsw.no> writes:

> Tim Rentsch wrote:
[...]
>> This technique of packing several small values into a constant
>> and shifting to get a specific result is useful in lots of
>> situations. For example, if we are writing a function to test
>> for primality, it can start with
>>
>> if( n < 64 ) return 0240404242024042424254 >>n &1;
>>
>> for a fast test on small integers.
>
> [...]
>
> Your 64-bit bitmap here is (inversely) related to the modulo 30
> packing I'm using in my sieve, where the 8 possible prime
> locations in each block of 30 numbers will fit in a byte.

I wrote one of these about 20 years, and recently was motivated
to write one again. I'm curious - have you tried running your
sieve for all primes less than a trillion (1e12)?

Terje Mathisen <terje.mathisen@tmsw.no> writes:

> Terje Mathisen wrote:
>
>> Since Tim is using a 4-way unroll, I believe he would be loading 4
>> values from a table of at least 7 entries to handle loop alignement.
>
> I took another look at this, and it looks pretty nice to use 4
> starting entries (a..d in increasing order), then a two-cycle
> recurrence:
>
> a = c+d
> b = c+2*d
> c = a+b
> d = a+2*b
>
> which every LEA-capable compiler around could convert into
>
> lea rax,[rcx+rdx]
> lea rbx,[rcx+2*rdx]
>
> lea rcx,[rax+rbx]
> lea rdx,[rax+2*rbx]
>
> The pairs of instructions can run at the same time, so this can run in
> two cycles as long as back-to-back LEA can run in subsequent cycles.
>
> To run even faster than this I seem to need actual multiplications,
> with small fibonacci numbers as the multipliers:
>
> f(n+1) = f(n)+f(n-1)
> f(n+2) = 2*f(n)+f(n-1)
> f(n+3) = 3*f(n)+2*f(n-1)
> f(n+4) = 5*f(n)+3*(f(n-1)
> ...
> f(n+8) = 34*(f(n)+21*f(n-1)

As much as I am tempted to post the actual code, let me start
with something more like a hint, which I fully expect you can
carry forward to do the whole deal.

(forward by 0) a b
(forward by 1) b a+b
(forward by 2) a+b a+2b
(forward by 3) a+2b 2a+3b
(forward by 4) 2a+3b 3a+5b
(forward by 5) 3a+5b 5a+8b
(forward by 6) 5a+8b 8a+13b

and so forth.

Terje Mathisen <terje.mathisen@tmsw.no> writes:

> Tim Rentsch wrote:

Just one short followup here for now...

>> By test, my O(log(n)) code runs in 25-30 nanoseconds for 128-bit
>
> That is impressive!

I took the O(log(n)) code and put it in python (with arbirary
precision integers). Running fib( 1000000 ) took well under
a second.

Exponentials are insanely slow, and O(log(n)) is insanely fast...

Tim Rentsch wrote:
> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>
>> Tim Rentsch wrote:
> [...]
>>> This technique of packing several small values into a constant
>>> and shifting to get a specific result is useful in lots of
>>> situations. For example, if we are writing a function to test
>>> for primality, it can start with
>>>
>>> if( n < 64 ) return 0240404242024042424254 >>n &1;
>>>
>>> for a fast test on small integers.
>>
>> [...]
>>
>> Your 64-bit bitmap here is (inversely) related to the modulo 30
>> packing I'm using in my sieve, where the 8 possible prime
>> locations in each block of 30 numbers will fit in a byte.
>
> I wrote one of these about 20 years, and recently was motivated
> to write one again. I'm curious - have you tried running your
> sieve for all primes less than a trillion (1e12)?
>
No, I wrote this back in the 32-bit days and even though it is blocked
(with a 128 kB block size corresponding to nearly 4M numbers sieved per
block afair), I never ran it _that_ far!

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

Terje Mathisen wrote:

>
> To run even faster than this I seem to need actual multiplications, with
> small fibonacci numbers as the multipliers:

> f(n+1) = f(n)+f(n-1)
> f(n+2) = 2*f(n)+f(n-1)
> f(n+3) = 3*f(n)+2*f(n-1)
> f(n+4) = 5*f(n)+3*(f(n-1)
> ....
> f(n+8) = 34*(f(n)+21*f(n-1)

Combining::

unsigned short table[18] = { 0, 1, 1, 2, 3, 5, 8, 13, 34, 55, 89,
144, 233, 377, 610, 987, 1597, 2584 };

fib( unsigned n )
{ i = n & 15;
a = table[i ];
b = table[i+1];
for( n -= i; n ; n -= 16 )
{
i = n & 15;
a = a*table[i ] + b*table[i+1];
b = a*table[i+1] + b*table[i+2];
}
return a;
}

Knocking off 16 fibs per loop.

MitchAlsup wrote:

> Terje Mathisen wrote:

>>
>> To run even faster than this I seem to need actual multiplications, with
>> small fibonacci numbers as the multipliers:

>> f(n+1) = f(n)+f(n-1)
>> f(n+2) = 2*f(n)+f(n-1)
>> f(n+3) = 3*f(n)+2*f(n-1)
>> f(n+4) = 5*f(n)+3*(f(n-1)
>> ....
>> f(n+8) = 34*(f(n)+21*f(n-1)

> Combining::

> unsigned short table[18] = { 0, 1, 1, 2, 3, 5, 8, 13, 34, 55, 89,
> 144, 233, 377, 610, 987, 1597, 2584 };

> fib( unsigned n )
> {
> i = n & 15;
> a = table[i ];
> b = table[i+1];
> for( n -= i; n ; n -= 16 )
> {
> i = n & 15;
t = a*table[i ] + b*table[i+1];
> b = a*table[i+1] + b*table[i+2];
a = t;
> }
> return a;
> }

> Knocking off 16 fibs per loop.

Tim Rentsch wrote:
> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>
>> Terje Mathisen wrote:
>>
>>> Since Tim is using a 4-way unroll, I believe he would be loading
>>> 4 values from a table of at least 7 entries to handle loop
>>> alignement.
>>
>> I took another look at this, and it looks pretty nice to use 4
>> starting entries (a..d in increasing order), then a two-cycle
>> recurrence:
>>
>> a = c+d b = c+2*d c = a+b d = a+2*b
>>
>> which every LEA-capable compiler around could convert into
>>
>> lea rax,[rcx+rdx] lea rbx,[rcx+2*rdx]
>>
>> lea rcx,[rax+rbx] lea rdx,[rax+2*rbx]
>>
>> The pairs of instructions can run at the same time, so this can run
>> in two cycles as long as back-to-back LEA can run in subsequent
>> cycles.
>>
>> To run even faster than this I seem to need actual
>> multiplications, with small fibonacci numbers as the multipliers:
>>
>> f(n+1) = f(n)+f(n-1) f(n+2) = 2*f(n)+f(n-1) f(n+3) =
>> 3*f(n)+2*f(n-1) f(n+4) = 5*f(n)+3*(f(n-1) ... f(n+8) =
>> 34*(f(n)+21*f(n-1)
>
> As much as I am tempted to post the actual code, let me start with
> something more like a hint, which I fully expect you can carry
> forward to do the whole deal.
>
> (forward by 0) a b (forward by 1) b a+b (forward
> by 2) a+b a+2b (forward by 3) a+2b 2a+3b (forward by
> 4) 2a+3b 3a+5b (forward by 5) 3a+5b 5a+8b (forward by 6)
> 5a+8b 8a+13b
>
> and so forth.

Isn't that exactly the same as what I posted just above, i.e. that
subsequent levels require pairs of multiplications by smaller fib numbers?

All these muls are of course independent, but since we only have a very
small number of multipliers in current CPUS, we cannot expect to do more
than a pair of them in each cycle (all with latency 3-5).

Doing it with adds instead means that we have to generate
a,2*a,3*a,5*a,8*a,13*a etc.

We can use LEA and shifts to handle a*(2,3,5,8) and b ditto, so that
gets us to the next 5 fib numbers all done in parallel and at latency 1
or 2:

All these can be done in parallel on a 5-wide CPU:
a5 = LEA(a,a*4);
b3 = LEA(b,b*2); b5 = LEA(b,b*4);
f(n) = a+b;
f(n+1) = LEA(b+a*2);

These next operations can run in the second cycle, so we've generated 5
new fib() numbers in two cycles:

f(n+2) = f(n)+f(n+1);
f(n+3) = a5+b3;
f(n+4) = LEA(b5+a*8);

In the same cycle we can generate a couple more aggregate multiples:

a13 = LEA(a5,a*8); b13=LEA(b5,b*8);

and then we get to spew out the 6th fib value in the third cycle

f(n+5) = LEA(a13,b*8);

but I get the feel that it is better to stay at 4 or 5 new numbers in
each step?

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

Terje Mathisen <terje.mathisen@tmsw.no> writes:

> Tim Rentsch wrote:
>
>> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>>
>>> Terje Mathisen wrote:
>>>
>>>> Since Tim is using a 4-way unroll, I believe he would be loading
>>>> 4 values from a table of at least 7 entries to handle loop
>>>> alignement.
>>>
>>> I took another look at this, and it looks pretty nice to use 4
>>> starting entries (a..d in increasing order), then a two-cycle
>>> recurrence:
>>>
>>> a = c+d b = c+2*d c = a+b d = a+2*b
>>>
>>> which every LEA-capable compiler around could convert into
>>>
>>> lea rax,[rcx+rdx] lea rbx,[rcx+2*rdx]
>>>
>>> lea rcx,[rax+rbx] lea rdx,[rax+2*rbx]
>>>
>>> The pairs of instructions can run at the same time, so this can run
>>> in two cycles as long as back-to-back LEA can run in subsequent
>>> cycles.
>>>
>>> To run even faster than this I seem to need actual
>>> multiplications, with small fibonacci numbers as the multipliers:
>>>
>>> f(n+1) = f(n)+f(n-1) f(n+2) = 2*f(n)+f(n-1) f(n+3) =
>>> 3*f(n)+2*f(n-1) f(n+4) = 5*f(n)+3*(f(n-1) ... f(n+8) =
>>> 34*(f(n)+21*f(n-1)
>>
>> As much as I am tempted to post the actual code, let me start with
>> something more like a hint, which I fully expect you can carry
>> forward to do the whole deal.
>>
>> (forward by 0) a b (forward by 1) b a+b (forward
>> by 2) a+b a+2b (forward by 3) a+2b 2a+3b (forward by
>> 4) 2a+3b 3a+5b (forward by 5) 3a+5b 5a+8b (forward by 6)
>> 5a+8b 8a+13b
>>
>> and so forth.
>
> Isn't that exactly the same as what I posted just above, i.e. that
> subsequent levels require pairs of multiplications by smaller fib
> numbers?

I didn't really follow what you were saying, but I'm pretty
sure the two aren't the same, because I keep only two values,
a and b, no matter how many levels of x-fold-speedup is done.

> All these muls are of course independent, but since we only have a
> very small number of multipliers in current CPUS, we cannot expect to
> do more than a pair of them in each cycle (all with latency 3-5).
>
> Doing it with adds instead means that we have to generate
> a,2*a,3*a,5*a,8*a,13*a etc.
>
> We can use LEA and shifts to handle a*(2,3,5,8) and b ditto, so that
> gets us to the next 5 fib numbers all done in parallel and at latency
> 1 or 2:
>
> All these can be done in parallel on a 5-wide CPU:
> a5 = LEA(a,a*4);
> b3 = LEA(b,b*2); b5 = LEA(b,b*4);
> f(n) = a+b;
> f(n+1) = LEA(b+a*2);
>
> These next operations can run in the second cycle, so we've generated
> 5 new fib() numbers in two cycles:
>
> f(n+2) = f(n)+f(n+1);
> f(n+3) = a5+b3;
> f(n+4) = LEA(b5+a*8);
>
> In the same cycle we can generate a couple more aggregate multiples:
>
> a13 = LEA(a5,a*8); b13=LEA(b5,b*8);
>
> and then we get to spew out the 6th fib value in the third cycle
>
> f(n+5) = LEA(a13,b*8);
>
> but I get the feel that it is better to stay at 4 or 5 new numbers in
> each step?

In the code that follows, b is the current fibonacci number, and a
is the previous fibonacci number.

U64
fib_example( unsigned n ){
U64 b = n&1, a = b^1;

if( n & 2 ) a += b, b += a;
if( n & 4 ) a += b, b += a, a += b, b += a;
if( n & 8 ){
U64 na = 13*a+21*b, nb = 21*a+34*b;
a = na, b = nb;
}

n >>= 4;
while( n-- ){
U64 na = 610*a + 987*b, nb = 987*a + 1597*b;
a = na, b = nb;
}

return b;
}

Note that fib(8) is 21, and fib(16) is 987.

It should be easy to see how to modify this code to extend it one
level and get a 32-fold speedup for the last stage, which does in
fact give a significant performance improvement.

Also I have no doubt that you can see ways to make the code run
even faster. I'm interested to see what you come up with, with
the stipulation that no lookup tables be used.

On 1/2/2024 2:08 AM, Tim Rentsch wrote:
> Daniel James <daniel@me.invalid> writes:
>
>> On 31/12/2023 09:12, Pancho wrote:
>>
>>> In fact, given that fib grows exponentially, so we only ever need to
>>> calculate it for small values of n, the only ?bad? solution is the
>>> naive double recursive one.
>>
>> Well, I might disagree.
>>
>> Given that we only ever need to calculate it for small values of n (as
>> you say, and if that's true) its efficiency doesn't matter. If you
>> need to compute it often enough that efficiency might become relevant
>> just put the results in a table and use that. Given that fib goes over
>> 10^12 after about 60 terms it needn't be a very big table.
>
> When the size of the table is more than ten times the size of
> the code, and the code is reasonably fast, using a table is a
> poor choice. Fibonacci functions don't have to be slow.
>
>> The naive doubly recursive solution has the huge merit of being
>> human-readable.
>
> But unacceptably slow.
>

Ironically, its slowness is its main use-case...

It puts function call/return through a workout, and also gives a simple
check value to detect if it has gone amiss.

There are many faster ways to calculate Fibonacci numbers, but not so
many that also happen to be a ready-made and fairly concise
microbenchmark...

Even if, yes, in premise a compiler could optimize it away.

But, say, if one wanted something less slow, say (untested):
long fasterfib(int n)
{ long c0, c1, c2;
if(n<2)
return(1);
c0=1;
c1=1;
c2=c0+c1;
while(n>2)
{
c0=c1;
c1=c2;
c2=c0+c1;
n--;
}
return(c2);
}

>> Unlike, say:
>>
>> double fib( int n )
>> {
>> return (pow( (1+sqrt(5.0))/2, n+1 )
>> - pow( (1-sqrt(5.0))/2, n+1 ))
>> / sqrt(5.0);
>> }
>>
>> ... which is great for fibs of big values, but not fast for small ones.
>
> Unacceptable for fibs of big values, because the values computed
> are wrong.

Tim Rentsch wrote:
> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>
>> Tim Rentsch wrote:
>>
>>> Terje Mathisen <terje.mathisen@tmsw.no> writes:
>>>
>>>> Terje Mathisen wrote:
>>>>
>>>>> Since Tim is using a 4-way unroll, I believe he would be loading
>>>>> 4 values from a table of at least 7 entries to handle loop
>>>>> alignement.
>>>>
>>>> I took another look at this, and it looks pretty nice to use 4
>>>> starting entries (a..d in increasing order), then a two-cycle
>>>> recurrence:
>>>>
>>>> a = c+d b = c+2*d c = a+b d = a+2*b
>>>>
>>>> which every LEA-capable compiler around could convert into
>>>>
>>>> lea rax,[rcx+rdx] lea rbx,[rcx+2*rdx]
>>>>
>>>> lea rcx,[rax+rbx] lea rdx,[rax+2*rbx]
>>>>
>>>> The pairs of instructions can run at the same time, so this can run
>>>> in two cycles as long as back-to-back LEA can run in subsequent
>>>> cycles.
>>>>
>>>> To run even faster than this I seem to need actual
>>>> multiplications, with small fibonacci numbers as the multipliers:
>>>>
>>>> f(n+1) = f(n)+f(n-1) f(n+2) = 2*f(n)+f(n-1) f(n+3) =
>>>> 3*f(n)+2*f(n-1) f(n+4) = 5*f(n)+3*(f(n-1) ... f(n+8) =
>>>> 34*(f(n)+21*f(n-1)
>>>
>>> As much as I am tempted to post the actual code, let me start with
>>> something more like a hint, which I fully expect you can carry
>>> forward to do the whole deal.
>>>
>>> (forward by 0) a b (forward by 1) b a+b (forward
>>> by 2) a+b a+2b (forward by 3) a+2b 2a+3b (forward by
>>> 4) 2a+3b 3a+5b (forward by 5) 3a+5b 5a+8b (forward by 6)
>>> 5a+8b 8a+13b
>>>
>>> and so forth.
>>
>> Isn't that exactly the same as what I posted just above, i.e. that
>> subsequent levels require pairs of multiplications by smaller fib
>> numbers?
>
> I didn't really follow what you were saying, but I'm pretty
> sure the two aren't the same, because I keep only two values,
> a and b, no matter how many levels of x-fold-speedup is done.
>
>> All these muls are of course independent, but since we only have a
>> very small number of multipliers in current CPUS, we cannot expect to
>> do more than a pair of them in each cycle (all with latency 3-5).
>>
>> Doing it with adds instead means that we have to generate
>> a,2*a,3*a,5*a,8*a,13*a etc.
>>
>> We can use LEA and shifts to handle a*(2,3,5,8) and b ditto, so that
>> gets us to the next 5 fib numbers all done in parallel and at latency
>> 1 or 2:
>>
>> All these can be done in parallel on a 5-wide CPU:
>> a5 = LEA(a,a*4);
>> b3 = LEA(b,b*2); b5 = LEA(b,b*4);
>> f(n) = a+b;
>> f(n+1) = LEA(b+a*2);
>>
>> These next operations can run in the second cycle, so we've generated
>> 5 new fib() numbers in two cycles:
>>
>> f(n+2) = f(n)+f(n+1);
>> f(n+3) = a5+b3;
>> f(n+4) = LEA(b5+a*8);
>>
>> In the same cycle we can generate a couple more aggregate multiples:
>>
>> a13 = LEA(a5,a*8); b13=LEA(b5,b*8);
>>
>> and then we get to spew out the 6th fib value in the third cycle
>>
>> f(n+5) = LEA(a13,b*8);
>>
>> but I get the feel that it is better to stay at 4 or 5 new numbers in
>> each step?
>
> In the code that follows, b is the current fibonacci number, and a
> is the previous fibonacci number.
>
> U64
> fib_example( unsigned n ){
> U64 b = n&1, a = b^1;
>
> if( n & 2 ) a += b, b += a;
> if( n & 4 ) a += b, b += a, a += b, b += a;
> if( n & 8 ){
> U64 na = 13*a+21*b, nb = 21*a+34*b;
> a = na, b = nb;
> }
>
> n >>= 4;
> while( n-- ){
> U64 na = 610*a + 987*b, nb = 987*a + 1597*b;
> a = na, b = nb;
> }
>
> return b;
> }
>
> Note that fib(8) is 21, and fib(16) is 987.

You are skipping intermediate numbers!

>
> It should be easy to see how to modify this code to extend it one
> level and get a 32-fold speedup for the last stage, which does in
> fact give a significant performance improvement.

I should have realized that you were in fact doing so since you showed
your log2(n) algorithm, but I was in fact assuming that you had to
generate all intermediate values as well.
>
> Also I have no doubt that you can see ways to make the code run
> even faster. I'm interested to see what you come up with, with
> the stipulation that no lookup tables be used.
>
The four MULs will be the limiting factor here, but the key idea is that
there will only ever be four of them, so this could go however far you
want to take it.

With four MULs you should get a new (a,b) pair every 6-7 clock cycles
(assuming two fully pipelined MUL units), this isn't that much faster
than code which generates 5 new fib values every two cycles, and then
you actually get to see all intermediate values as well.

Still very nice idea. :-)

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

mitchalsup@aol.com (MitchAlsup) writes:

> MitchAlsup wrote:
>
>> Terje Mathisen wrote:
>
>>>
>>> To run even faster than this I seem to need actual multiplications,
>>> with small fibonacci numbers as the multipliers:
>>>
>>> f(n+1) = f(n)+f(n-1)
>>> f(n+2) = 2*f(n)+f(n-1)
>>> f(n+3) = 3*f(n)+2*f(n-1)
>>> f(n+4) = 5*f(n)+3*(f(n-1)
>>> ....
>>> f(n+8) = 34*(f(n)+21*f(n-1)
>>
>> Combining::
>>
>> unsigned short table[18] = { 0, 1, 1, 2, 3, 5, 8, 13, 34, 55, 89,
>> 144, 233, 377, 610, 987, 1597, 2584 };
>>
>> fib( unsigned n )
>> {
>> i = n & 15;
>> a = table[i ];
>> b = table[i+1];
>> for( n -= i; n ; n -= 16 )
>> {
>> i = n & 15;
>
> t = a*table[i ] + b*table[i+1];
>
>> b = a*table[i+1] + b*table[i+2];
>
> a = t;

After n -= i, the lower order four bits of n are zero, which means
the coefficients being used are 0, 1, and 1. Also the table is
missing the value 21.

U64
fib( unsigned n ){
static unsigned short table[17] = {
1, 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610,
};
U64 a = table[n&15], b = table[(n&15)+1];

for( n >>= 4; n--; ){
U64 na = 610*a + 987*b, nb = 987*a + 1597*b;
a = na, b = nb;
}

return b;
}

BGB <cr88192@gmail.com> writes:

> On 1/2/2024 2:08 AM, Tim Rentsch wrote:
>
>> Daniel James <daniel@me.invalid> writes:
>>
>>> On 31/12/2023 09:12, Pancho wrote:
>>>
>>>> In fact, given that fib grows exponentially, so we only ever need to
>>>> calculate it for small values of n, the only ?bad? solution is the
>>>> naive double recursive one.
>>>
>>> Well, I might disagree.
>>>
>>> Given that we only ever need to calculate it for small values of n (as
>>> you say, and if that's true) its efficiency doesn't matter. If you
>>> need to compute it often enough that efficiency might become relevant
>>> just put the results in a table and use that. Given that fib goes over
>>> 10^12 after about 60 terms it needn't be a very big table.
>>
>> When the size of the table is more than ten times the size of
>> the code, and the code is reasonably fast, using a table is a
>> poor choice. Fibonacci functions don't have to be slow.
>>
>>> The naive doubly recursive solution has the huge merit of being
>>> human-readable.
>>
>> But unacceptably slow.

[..]

> But, say, if one wanted something less slow, say (untested):
> long fasterfib(int n)
> {
> long c0, c1, c2;
> if(n<2)
> return(1);
> c0=1;
> c1=1;
> c2=c0+c1;
> while(n>2)
> {
> c0=c1;
> c1=c2;
> c2=c0+c1;
> n--;
> }
> return(c2);
> }

The Fibonacci numbers start at 0: fib(0) == 0, fib(1) == 1, ...
The code above produces something else.

What you're showing is a standard iterative version, much like
the one posted by Thomas Koenig early on.

Terje Mathisen <terje.mathisen@tmsw.no> writes:

> Tim Rentsch wrote:

[...]

>> In the code that follows, b is the current fibonacci number, and a
>> is the previous fibonacci number.
>>
>> U64
>> fib_example( unsigned n ){
>> U64 b = n&1, a = b^1;
>>
>> if( n & 2 ) a += b, b += a;
>> if( n & 4 ) a += b, b += a, a += b, b += a;
>> if( n & 8 ){
>> U64 na = 13*a+21*b, nb = 21*a+34*b;
>> a = na, b = nb;
>> }
>>
>> n >>= 4;
>> while( n-- ){
>> U64 na = 610*a + 987*b, nb = 987*a + 1597*b;
>> a = na, b = nb;
>> }
>>
>> return b;
>> }
>>
>> Note that fib(8) is 21, and fib(16) is 987.
>
> You are skipping intermediate numbers!

Oh yes. It never occurred to me to do anything but.

[...]

> Still very nice idea. :-)

Thank you. It seem obvious in retrospect, but that property
holds for lots of ideas that are not obvious ahead of time.

I have an interesting footnote. My first implementation of this
idea was written using a functional recursive style. I was doing
performance measurements using a relatively low optimization setting
(-Os). Just this morning I decided to turn up the level to -O3, and
was surprised to see that the functional version ran faster than
everything else, 25% faster than the closest imperative version.
I've seen this result before in other cases - code written in a
functional style actually runs faster than code written in an
imperative style.

Here is the code:

typedef unsigned long long U64;

static U64 ff1( U64, U64, unsigned );
static U64 ff2( U64, U64, unsigned );
static U64 ff3( U64, U64, unsigned );
static U64 ff4( U64, U64, unsigned );

U64
functional_fib( unsigned n ){
return
n < 8 ? 0xd8532110u >>(n*4) &15 :
n & 1 ? ff1( 0, 1, n/2 ) :
/*******/ ff1( 1, 0, n/2 );
}

U64
ff1( U64 a, U64 b, unsigned k ){
return k & 1 ? ff2( a+b, a+2*b, k/2 ) : ff2( a, b, k/2 );
}

U64
ff2( U64 a, U64 b, unsigned k ){
return k & 1 ? ff3( 2*a+3*b, 3*a+5*b, k/2 ) : ff3( a, b, k/2 );
}

U64
ff3( U64 a, U64 b, unsigned k ){
return k & 1 ? ff4( 13*a+21*b, 21*a+34*b, k/2 ) : ff4( a, b, k/2 );
}

U64
ff4( U64 a, U64 b, unsigned k ){
return k ? ff4( 610*a+987*b, 987*a+1597*b, k-1 ) : b;
}

On 1/2/24 2:10 PM, MitchAlsup wrote:
> Thomas Koenig wrote:
>
>> Compilers often do not do a good job of replacing recursion
>> with iteration (is there a verb for that? Maybe it can be called
>> iterationizing?), so it is something to consider if an algorithm
>> is at all time critical or stack space is problematic.
>
>
> Tail recursion elimination.

I thought compilers generally were good at tail recursion
elimination but did not attempt "flattening" something like
flood-fill (which is not "tail recursion" but could avoid
some of the function call overhead). (I.e., I vaguely
remember over a decade ago noticing that this case was not
optimized — for the samegnome game which has a small board
so recursive depth and call overhead was not really significant.)

I also vaguely recall reading at some point that even
recursion that was "easy" to convert to tail recursion was
not handled by gcc (also many years ago).