[PATCH] Fibonacci/Lucas tweaks

[marco.bodrato at tutanota.com]

Ciao David,

13 feb 2026, 05:36 da sparks05 at proton.me:

> Here's the results of my poking at the Lucas number evaluation code.
>
> ----8<----
> Patch 1/3: mpz_lucnum_ui: Use a simpler Lucas formula.
>  
> The code previously used
> L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*{-1)^k
>  = 5*F[k-1]*F[k+2]          - 4*{-1)^k
> But it's much simpler to ask for the next Fibonacci pair and use
> L[2k+1] = 5*F[k]*F[k+1] + (-1)^k 
>
You are right!
It is known that
L[k+m] = F[k]*F[m]+ L[m-k](-1)^k

Using m=k+1, and L[1]=1 as you propose is much simpler and straightforward
than using F[k-1] and F[k+2]=2F[k]+F[k-1], and L[k+2-(k-1)]=L[3]=4.

> The diff is larger that expected because this also swaps the
> order of tests in the loop, so the test for the FIB_TABLE comes first.
>
Why did you swap that test?

The current code:
 - check if n is smaller than the table limit
 - then loop and in the loop
   + checks if n is odd;
   + otherwise divides by 2, and
   + checks if the divided n is smaller than the table limit.

Checking at the beginning should not be worth it, because we just checked before the loop, didn't we? When we start the loop we already know that n does not fit as a table index, we can wait after the division before checking again.
 

> The first hunk shrinks the temporary buffer size slightly, saving
> one limb much of the time.
>
Seems like a good idea.

> +	  /* L[n] = F[n] + 2F[n-1] */
> +	  lp[0] = FIB_TABLE (n) + 2 * FIB_TABLE ((int) n - 1);
>
I wonder if it's wort adding a condition:
  if (FIB_TABLE_LUCNUM_LIMIT < FIB_TABLE_LIMIT)
    /* L[n] = F[n+1] + F[n-1] */
    lp[0] = FIB_TABLE ((int) n + 1) + FIB_TABLE ((int) n - 1);
  else
    /* L[n] = F[n] + 2F[n-1] */
    lp[0] = FIB_TABLE (n) + 2 * FIB_TABLE ((int) n - 1);

This is not a branch at runtime, because the condition can be checked at compile time.
But, is it better to read two adjacent values from memory or to avoid a not expensive shift?
I'm not sure about the answer.

> ----8<----
> Patch 2/3: MPN_FIB2_SIZE: Prepare for tighter approximation.
>
Let's look at this.

> diff --git a/mpz/lucnum2_ui.c b/mpz/lucnum2_ui.c
> index 3cc3d7b75..554aba0b0 100644
> --- a/mpz/lucnum2_ui.c
> +++ b/mpz/lucnum2_ui.c
> @@ -61,7 +61,7 @@ mpz_lucnum2_ui (mpz_ptr ln, mpz_ptr lnsub1, unsigned long n)
>  }
>  
>  TMP_MARK;
> -  size = MPN_FIB2_SIZE (n);
> +  size = MPN_FIB2_SIZE (n+2);
>  f1p = TMP_ALLOC_LIMBS (size);
>
You are right!

>  /* check MPN_FIB2_SIZE seems right, compared to actual size and
> -         compared to our float formula */
> -      if (MPN_FIB2_SIZE (n) < MPZ_FIB_SIZE_FLOAT (n))
> +	 compared to our float formula.  MPN_FIB2_SIZE has a known problem
> +	 when n == 1 (it returns 0 rather than 1) which doesn't affect
> +	 any real caller, so we accomodate it in the test. */
>
Does MPN_FIB2_SIZE really have that problem?
The 1+ at the beginning should address the issue, doesn't it?

> ----8<----
> Patch 3/3:    MPN_FIB2_SIZE: Tighten approximation.
>
Ok, let's look at it.
> -#define MPN_FIB2_SIZE(n) \
> -  ((mp_size_t) ((n) / 32 * 23 / GMP_NUMB_BITS) + 4)
>
The old code needed: a couple of shifts, a multiplication, an addition.
> +#define MPN_FIB2_SIZE(n) (1 + \
> +   (mp_size_t)( (n)    / 65536 * MPN_FIB2_MULT) + \
> +   (mp_size_t)(((n)|1) % 65536 * MPN_FIB2_MULT + 65536-MPN_FIB2_SUB) / 65536)
>
The proposed one: a couple of shifts, an and mask, 3 additions.
How much does the MPN_FIB2_SUB term change the result?
Should we replace this estimate with a single umul_ppmm (_result, _dummy, n, _const) like we do in macros like DIGITS_IN_BASE_FROM_BITS?

The constant could be generated by gen-fib.c

Ĝis,
m
-------------- next part --------------
A non-text attachment was scrubbed...
Name: gen-fib_limbsPerIter.patch
Type: text/x-patch
Size: 1663 bytes
Desc: not available
URL: <https://gmplib.org/list-archives/gmp-devel/attachments/20260217/5b467349/attachment-0001.bin>