RosettaCodeData/Task/Fibonacci-sequence/Pascal/fibonacci-sequence-1.pas

Program FastFibonacciGMP;
{$mode objfpc}{$H+}{$J-}

// ============================================================================
// Fast Fibonacci Calculator with Multiple Algorithms and GMP
// ----------------------------------------------------------------------------
// 2025.12.05
// Author: jpd
// AI Assistant: DeepSeek
// Description: Multiple Fibonacci algorithms with arbitrary precision
//              AI translated from my own 2016 python app ( the hybrid part )
//
// LICENSE & ATTRIBUTION:
// ----------------------------------------------------------------------------
// 1. This Pascal source file (`FastFibonacciGMP.pas`) is released to the
//    PUBLIC DOMAIN for Rosetta Code contribution. Author waives all copyright.
//
// 2. This program uses the Free Pascal `GMP` unit, which is a binding to the
//    GNU MP Library (GMP). The `GMP` unit is part of Free Pascal's RTL-extra
//    package, licensed under the GNU Lesser General Public License (LGPL)
//    with a static linking exception.
//
// 3. The underlying GNU MP Library (libgmp) itself is licensed under:
//    - GNU Lesser General Public License version 3 or later (LGPLv3+), OR
//    - GNU General Public License version 2 or later (GPLv2+).
//
// 4. When compiled, this program links to the external GMP C library.
//    For license compliance, users must have access to the GMP library
//    source code (available at https://gmplib.org/).
// ============================================================================
// MATHEMATICAL FOUNDATION ATTRIBUTION:
// ----------------------------------------------------------------------------
// 1. Fast Doubling Formulas (Fibonacci Squaring):
//    F(2k) = F(k) × [2×F(k+1) - F(k)]
//    F(2k+1) = F(k+1)² + F(k)²
//    Derived from Cassini's identity (1680) and Catalan's identity (1879).
//    Modern algorithmic presentation: Dijkstra (1978), Cohn (1963).
//
// 2. Fibonacci Tripling Formulas:
//    F(3k) = 5×F(k)³ + 3×(-1)ᵏ×F(k)
//    F(3k+1) = F(k+1)³ + 3×F(k+1)×F(k)² - F(k)³
//    Derived from Binet's formula (1843).
//    Algorithmic optimization: Various number theory sources.
//
// 3. Hybrid Decomposition Algorithm:
//    Recursive decomposition using factors 2 and 3 based on divisibility.
//    Original implementation concept: jpd (2016 Python implementation).
//    Ported to Pascal with algorithmic corrections: DeepSeek AI (2025).
//
// NOTE: These mathematical identities are in the public domain.
//       This specific implementation is original work.
// ----------------------------------------------------------------------------
Uses
    SysUtils,
    DateUtils,
    Gmp;

// ----------------------------------------------------------------------------
// Matrix type for 2x2 integer matrix using GMP arbitrary precision
// ----------------------------------------------------------------------------
Type
    TMatrix = Record
        a, b, c, d: mpz_t;
    End;

    // Enum for algorithm selection
    TAlgorithm = (algMatrix, algFastDoubling, algIterative, algHybrid);

// ----------------------------------------------------------------------------
// Initialize matrix with four integer values
// ----------------------------------------------------------------------------
Procedure InitMatrix(Var m: TMatrix; a, b, c, d: Integer);
Begin
    mpz_init_set_si(m.a, a);
    mpz_init_set_si(m.b, b);
    mpz_init_set_si(m.c, c);
    mpz_init_set_si(m.d, d);
End;

// ----------------------------------------------------------------------------
// Clear GMP memory for matrix elements
// ----------------------------------------------------------------------------
Procedure ClearMatrix(Var m: TMatrix);
Begin
    mpz_clear(m.a);
    mpz_clear(m.b);
    mpz_clear(m.c);
    mpz_clear(m.d);
End;

// ----------------------------------------------------------------------------
// Multiply two matrices: R = A × B
// Uses GMP for arbitrary precision arithmetic
// ----------------------------------------------------------------------------
Procedure MatrixMul(Var R, A, B: TMatrix);
Var
    temp1, temp2, temp3, temp4: mpz_t;
Begin
    // Initialize temporaries
    mpz_init(temp1);
    mpz_init(temp2);
    mpz_init(temp3);
    mpz_init(temp4);

    // temp1 = A.a × B.a + A.b × B.c
    mpz_mul(temp1, A.a, B.a);
    mpz_mul(temp2, A.b, B.c);
    mpz_add(temp1, temp1, temp2);

    // temp2 = A.a × B.b + A.b × B.d
    mpz_mul(temp2, A.a, B.b);
    mpz_mul(temp3, A.b, B.d);
    mpz_add(temp2, temp2, temp3);

    // temp3 = A.c × B.a + A.d × B.c
    mpz_mul(temp3, A.c, B.a);
    mpz_mul(temp4, A.d, B.c);
    mpz_add(temp3, temp3, temp4);

    // temp4 = A.c × B.b + A.d × B.d
    mpz_mul(temp4, A.c, B.b);
    mpz_mul(R.d, A.d, B.d);  // Reuse R.d as temp
    mpz_add(temp4, temp4, R.d);

    // Store results
    mpz_set(R.a, temp1);
    mpz_set(R.b, temp2);
    mpz_set(R.c, temp3);
    mpz_set(R.d, temp4);

    // Cleanup temporaries
    mpz_clear(temp1);
    mpz_clear(temp2);
    mpz_clear(temp3);
    mpz_clear(temp4);
End;

// ----------------------------------------------------------------------------
// Compute Fibonacci number F(n) using matrix exponentiation
// Supports both positive and negative indices using 64-bit integers
// ----------------------------------------------------------------------------
Procedure FibonacciMatrix(n: Int64; Var result: mpz_t);
Var
    Base, ResultMat, TempMat: TMatrix;
    mask: QWord;
    abs_n: QWord;
    startTime, endTime: TDateTime;
Begin
    startTime := Now;

    // Initialize matrices
    If n >= 0 Then
        InitMatrix(Base, 0, 1, 1, 1)      // [0 1; 1 1] for positive n
    Else
        InitMatrix(Base, -1, 1, 1, 0);    // [-1 1; 1 0] for negative n

    // Identity matrix for exponentiation
    InitMatrix(ResultMat, 1, 0, 0, 1);
    InitMatrix(TempMat, 0, 0, 0, 0);

    // Use absolute value for exponentiation
    If n >= 0 Then
        abs_n := QWord(n)
    Else
        abs_n := QWord(-n);

    mask := 1;

    // Fast exponentiation by squaring: O(log n) complexity
    While mask <= abs_n Do
    Begin
        If (abs_n And mask) <> 0 Then
        Begin
            MatrixMul(TempMat, ResultMat, Base);
            mpz_set(ResultMat.a, TempMat.a);
            mpz_set(ResultMat.b, TempMat.b);
            mpz_set(ResultMat.c, TempMat.c);
            mpz_set(ResultMat.d, TempMat.d);
        End;

        // Square the base matrix
        MatrixMul(TempMat, Base, Base);
        mpz_set(Base.a, TempMat.a);
        mpz_set(Base.b, TempMat.b);
        mpz_set(Base.c, TempMat.c);
        mpz_set(Base.d, TempMat.d);

        mask := mask Shl 1;
    End;

    // Result is in the (0,1) position of the matrix
    mpz_set(result, ResultMat.b);

    // Cleanup
    ClearMatrix(Base);
    ClearMatrix(ResultMat);
    ClearMatrix(TempMat);

    endTime := Now;
    WriteLn('Matrix method time: ', MilliSecondsBetween(endTime, startTime), ' ms');
End;

// ----------------------------------------------------------------------------
// Iterative version of Fast Doubling (avoids recursion depth issues)*
//
// Process bits of n from most significant to least significant
// ----------------------------------------------------------------------------
Procedure FastDoublingIterative(n: QWord; Var result: mpz_t);
Var
    a, b, c, d, temp1, temp2: mpz_t;
    mask: QWord;
Begin
    // Initialize GMP variables
    mpz_init_set_ui(a, 0);  // F(0) = 0
    mpz_init_set_ui(b, 1);  // F(1) = 1
    mpz_init(c);
    mpz_init(d);
    mpz_init(temp1);
    mpz_init(temp2);

    // Find highest set bit
    mask := 1;
    While mask <= n Do
        mask := mask Shl 1;
    mask := mask Shr 1;  // Now mask is the highest set bit

    // Process bits from most significant to least significant
    While mask > 0 Do
    Begin
        // c = F(2k) = F(k) * [2 * F(k+1) - F(k)]
        // where k = current value of a,b (representing F(k), F(k+1))

        // Calculate 2*b - a
        mpz_mul_ui(temp1, b, 2);    // 2 * F(k+1)
        mpz_sub(temp1, temp1, a);   // 2 * F(k+1) - F(k)

        // c = a * (2*b - a)
        mpz_mul(c, a, temp1);

        // d = F(k+1)² + F(k)²
        mpz_mul(temp1, b, b);       // b²
        mpz_mul(temp2, a, a);       // a²
        mpz_add(d, temp1, temp2);   // a² + b²

        // Update a, b based on current bit of n
        If (n And mask) <> 0 Then
        Begin
            // Current bit is 1: (a,b) = (d, c+d)
            // This corresponds to going from (F(k), F(k+1)) to (F(2k+1), F(2k+2))
            mpz_add(temp1, c, d);   // c + d = F(2k+1) + F(2k) = F(2k+2)
            mpz_set(a, d);          // a = F(2k+1)
            mpz_set(b, temp1);      // b = F(2k+2)
        End
        Else
        Begin
            // Current bit is 0: (a,b) = (c, d)
            // This corresponds to going from (F(k), F(k+1)) to (F(2k), F(2k+1))
            mpz_set(a, c);          // a = F(2k)
            mpz_set(b, d);          // b = F(2k+1)
        End;

        mask := mask Shr 1;
    End;

    // Result is in a
    mpz_set(result, a);

    // Cleanup
    mpz_clear(a);
    mpz_clear(b);
    mpz_clear(c);
    mpz_clear(d);
    mpz_clear(temp1);
    mpz_clear(temp2);
End;

// ----------------------------------------------------------------------------
// Wrapper for Fast Doubling with timing
// ----------------------------------------------------------------------------
Procedure FibonacciFastDoubling(n: Int64; Var result: mpz_t);
Var
    startTime, endTime: TDateTime;
    fn: mpz_t;
    abs_n: QWord;
Begin
    startTime := Now;

    // Handle negative n using F(-n) = (-1)^(n+1) * F(n)
    If n >= 0 Then
        abs_n := QWord(n)
    Else
        abs_n := QWord(-n);

    mpz_init(fn);
    FastDoublingIterative(abs_n, fn);

    If n < 0 Then
    Begin
        // Apply sign for negative indices
        If (abs_n Mod 2) = 0 Then
            mpz_neg(fn, fn);
    End;

    mpz_set(result, fn);
    mpz_clear(fn);

    endTime := Now;
    WriteLn('Fast doubling method time: ', MilliSecondsBetween(endTime, startTime), ' ms');
End;

// ----------------------------------------------------------------------------
// Hybrid Doubling/Tripling Method (Fastest and most efficient)
// Based on Python matfib.1.1.py algorithm with corrected logic*
// ----------------------------------------------------------------------------
Procedure FibonacciHybrid(n: Int64; Var result: mpz_t);
Var
    startTime, endTime: TDateTime;
    FDict: array of record
        key: QWord;
        value: mpz_t;
    end;
    DivMap: array of record
        key: QWord;
        divisor: QWord;
    end;
    i, j, idx: Integer;
    a, divBy, nextNum: QWord;
    temp1, temp2, temp3: mpz_t;
    sign: Boolean;
    abs_n: QWord;
    target: QWord;
    sortedKeys: array of QWord;

    // Find index of key in FDict array
    Function FindF(key: QWord): Integer;
    Var
        i: Integer;
    Begin
        For i := 0 To Length(FDict) - 1 Do
            If FDict[i].key = key Then
            Begin
                Result := i;
                Exit;
            End;
        Result := -1;
    End;

    // Add new Fibonacci number to FDict array
    Procedure AddF(key: QWord);
    Begin
        SetLength(FDict, Length(FDict) + 1);
        FDict[Length(FDict) - 1].key := key;
        mpz_init(FDict[Length(FDict) - 1].value);
    End;

    // Helper: Compute F(n+1) if not already in dictionary*
    Procedure EnsureFPlusOne(n: QWord);
    Var
        idx1, idx2: Integer;
    Begin
        If FindF(n + 1) >= 0 Then
            Exit;

        // Compute F(n+1) = F(n) + F(n-1)
        idx1 := FindF(n);
        idx2 := FindF(n - 1);

        If (idx1 < 0) or (idx2 < 0) Then
        Begin
            WriteLn('Error: Cannot compute F(', n+1, ')');
            Halt(1);
        End;

        AddF(n + 1);
        mpz_add(FDict[Length(FDict) - 1].value, FDict[idx1].value, FDict[idx2].value);
    End;

Begin
    startTime := Now;

    // Initialize GMP temporaries
    mpz_init(temp1);
    mpz_init(temp2);
    mpz_init(temp3);

    // Handle negative indices
    If n >= 0 Then
        abs_n := QWord(n)
    Else
        abs_n := QWord(-n);

    target := abs_n;
    sign := (n < 0) And ((abs_n Mod 2) = 0);

    // Initialize base Fibonacci numbers (0-8)
    SetLength(FDict, 9);
    For i := 0 To 8 Do
    Begin
        FDict[i].key := QWord(i);
        mpz_init(FDict[i].value);
    End;

    mpz_set_ui(FDict[0].value, 0);   // F(0) = 0
    mpz_set_ui(FDict[1].value, 1);   // F(1) = 1
    mpz_set_ui(FDict[2].value, 1);   // F(2) = 1
    mpz_set_ui(FDict[3].value, 2);   // F(3) = 2
    mpz_set_ui(FDict[4].value, 3);   // F(4) = 3
    mpz_set_ui(FDict[5].value, 5);   // F(5) = 5
    mpz_set_ui(FDict[6].value, 8);   // F(6) = 8
    mpz_set_ui(FDict[7].value, 13);  // F(7) = 13
    mpz_set_ui(FDict[8].value, 21);  // F(8) = 21

    // If n is small, return directly
    If abs_n <= 8 Then
    Begin
        idx := FindF(abs_n);
        If idx >= 0 Then
            mpz_set(result, FDict[idx].value);

        If sign Then
            mpz_neg(result, result);

        // Cleanup
        For i := 0 To 8 Do
            mpz_clear(FDict[i].value);
        mpz_clear(temp1);
        mpz_clear(temp2);
        mpz_clear(temp3);

        endTime := Now;
        WriteLn('Hybrid method time: ', MilliSecondsBetween(endTime, startTime), ' ms');
        Exit;
    End;

    // Build division path
    SetLength(DivMap, 0);
    a := abs_n;

    // Python: while a > 1:
    While a > 1 Do
    Begin
        // Python: div = 3 if a % 3.0 == 0.0 else 2
        If (a Mod 3) = 0 Then
            divBy := 3
        Else
            divBy := 2;

        a := a Div divBy;

        // Python: _[a] = div
        SetLength(DivMap, Length(DivMap) + 1);
        DivMap[Length(DivMap) - 1].key := a;
        DivMap[Length(DivMap) - 1].divisor := divBy;
    End;

    // Python: _[a] = div (final assignment)
    SetLength(DivMap, Length(DivMap) + 1);
    DivMap[Length(DivMap) - 1].key := a;
    DivMap[Length(DivMap) - 1].divisor := divBy;

    // Sort keys in increasing order
    SetLength(sortedKeys, Length(DivMap));
    For i := 0 To Length(DivMap) - 1 Do
        sortedKeys[i] := DivMap[i].key;

    // Simple bubble sort
    For i := 0 To Length(sortedKeys) - 2 Do
        For j := i + 1 To Length(sortedKeys) - 1 Do
            If sortedKeys[i] > sortedKeys[j] Then
            Begin
                a := sortedKeys[i];
                sortedKeys[i] := sortedKeys[j];
                sortedKeys[j] := a;
            End;

    Write('Operations: ');

    // Process each key in sorted order
    For i := 0 To Length(sortedKeys) - 1 Do
    Begin
        a := sortedKeys[i];

        // Find the divisor for this key
        divBy := 0;
        For j := 0 To Length(DivMap) - 1 Do
            If DivMap[j].key = a Then
            Begin
                divBy := DivMap[j].divisor;
                Break;
            End;

        If divBy = 0 Then
            Continue;

        nextNum := a * divBy;  // Python: i = _[n] * n

        If divBy = 2 Then
        Begin
            Write('² ');

            // Need F(a-1) and F(a)
            j := FindF(a - 1);
            idx := FindF(a);

            If (j < 0) or (idx < 0) Then
            Begin
                WriteLn('Error: Missing base Fibonacci numbers for doubling at a=', a);
                Halt(1);
            End;

            // F(i-1) = F(a-1)² + F(a)²
            AddF(nextNum - 1);
            mpz_mul(temp1, FDict[j].value, FDict[j].value);  // F(a-1)²
            mpz_mul(temp2, FDict[idx].value, FDict[idx].value);  // F(a)²
            mpz_add(FDict[Length(FDict) - 1].value, temp1, temp2);

            // Check if we found target
            If (nextNum - 1) = target Then
                Break;

            // F(i) = (2 * F(a-1) + F(a)) * F(a)
            AddF(nextNum);
            mpz_mul_ui(temp1, FDict[j].value, 2);  // 2 * F(a-1)
            mpz_add(temp1, temp1, FDict[idx].value);  // 2*F(a-1) + F(a)
            mpz_mul(FDict[Length(FDict) - 1].value, temp1, FDict[idx].value);

            // Check if we found target
            If nextNum = target Then
                Break;

            // CRITICAL FIX: Always compute F(i+1) and F(i+2)**
            // F(i+1) = F(i-1) + F(i)
            AddF(nextNum + 1);
            j := FindF(nextNum - 1);
            idx := FindF(nextNum);
            mpz_add(FDict[Length(FDict) - 1].value, FDict[j].value, FDict[idx].value);

            // Check if we found target
            If (nextNum + 1) = target Then
                Break;

            // F(i+2) = F(i) + F(i+1)
            AddF(nextNum + 2);
            j := FindF(nextNum);
            idx := FindF(nextNum + 1);
            mpz_add(FDict[Length(FDict) - 1].value, FDict[j].value, FDict[idx].value);

            // Check if we found target
            If (nextNum + 2) = target Then
                Break;
        End
        Else If divBy = 3 Then
        Begin
            Write('³ ');

            // Need F(a) and ensure F(a+1) is available**
            idx := FindF(a);
            If idx < 0 Then
            Begin
                WriteLn('Error: Missing base Fibonacci number for tripling at a=', a);
                Halt(1);
            End;

            // Ensure F(a+1) is available**
            EnsureFPlusOne(a);
            j := FindF(a + 1);

            // F(i) = 5 * F(a)³ + 3 * (-1)^a * F(a)
            AddF(nextNum);

            // Compute F(a)² and F(a)³
            mpz_mul(temp1, FDict[idx].value, FDict[idx].value);  // F(a)²
            mpz_mul(temp2, temp1, FDict[idx].value);  // F(a)³

            // 5 * F(a)³
            mpz_mul_ui(temp3, temp2, 5);

            // 3 * (-1)^a * F(a)
            If (a Mod 2) = 0 Then
            Begin
                // (-1)^a = 1 for even a
                mpz_mul_ui(temp1, FDict[idx].value, 3);
                mpz_add(FDict[Length(FDict) - 1].value, temp3, temp1);
            End
            Else
            Begin
                // (-1)^a = -1 for odd a
                mpz_mul_ui(temp1, FDict[idx].value, 3);
                mpz_neg(temp1, temp1);
                mpz_add(FDict[Length(FDict) - 1].value, temp3, temp1);
            End;

            // Check if we found target
            If nextNum = target Then
                Break;

            // CRITICAL FIX: Always compute the window**
            // F(i+1) = F(a+1)³ + 3 * F(a+1) * F(a)² - F(a)³
            AddF(nextNum + 1);

            // F(a+1)³
            mpz_mul(temp1, FDict[j].value, FDict[j].value);
            mpz_mul(temp1, temp1, FDict[j].value);

            // 3 * F(a+1) * F(a)²
            mpz_mul(temp2, FDict[idx].value, FDict[idx].value);  // F(a)²
            mpz_mul_ui(temp3, FDict[j].value, 3);  // 3 * F(a+1)
            mpz_mul(temp3, temp3, temp2);  // 3 * F(a+1) * F(a)²

            // F(a)³
            mpz_mul(temp2, temp2, FDict[idx].value);  // F(a)³

            // F(i+1) = F(a+1)³ + 3*F(a+1)*F(a)² - F(a)³
            mpz_add(FDict[Length(FDict) - 1].value, temp1, temp3);
            mpz_sub(FDict[Length(FDict) - 1].value, FDict[Length(FDict) - 1].value, temp2);

            // Check if we found target
            If (nextNum + 1) = target Then
                Break;

            // F(i-1) = F(i+1) - F(i)
            AddF(nextNum - 1);
            j := FindF(nextNum + 1);
            idx := FindF(nextNum);
            mpz_sub(FDict[Length(FDict) - 1].value, FDict[j].value, FDict[idx].value);

            // Check if we found target
            If (nextNum - 1) = target Then
                Break;

            // F(i+2) = F(i) + F(i+1)
            AddF(nextNum + 2);
            j := FindF(nextNum);
            idx := FindF(nextNum + 1);
            mpz_add(FDict[Length(FDict) - 1].value, FDict[j].value, FDict[idx].value);

            // Check if we found target
            If (nextNum + 2) = target Then
                Break;

            // F(i+3) = F(i+1) + F(i+2)
            AddF(nextNum + 3);
            j := FindF(nextNum + 1);
            idx := FindF(nextNum + 2);
            mpz_add(FDict[Length(FDict) - 1].value, FDict[j].value, FDict[idx].value);

            // Check if we found target
            If (nextNum + 3) = target Then
                Break;
        End;
    End;

    WriteLn;

    // Get the result
    idx := FindF(target);
    If idx < 0 Then
    Begin
        WriteLn('Error: Could not compute F(', n, ')');
        Halt(1);
    End;

    mpz_set(result, FDict[idx].value);

    // Apply sign for negative indices
    If sign Then
        mpz_neg(result, result);

    // Cleanup
    For i := 0 To Length(FDict) - 1 Do
        mpz_clear(FDict[i].value);

    mpz_clear(temp1);
    mpz_clear(temp2);
    mpz_clear(temp3);

    endTime := Now;
    WriteLn('Hybrid method time: ', MilliSecondsBetween(endTime, startTime), ' ms');
End;

// ----------------------------------------------------------------------------
// Display Fibonacci number with first/last digits and total count
// ----------------------------------------------------------------------------
Procedure DisplayLargeFibonacci(n: Int64; Var fib: mpz_t);
Var
    s: AnsiString;
    len: Integer;
    first20, last20: AnsiString;
Begin
    s := mpz_get_str(Nil, 10, fib);
    len := Length(s);

    WriteLn('F(', n, ') has ', len, ' digits');

    If len <= 40 Then
        WriteLn('F(', n, ') = ', s)
    Else
    Begin
        first20 := Copy(s, 1, 20);
        last20 := Copy(s, len - 19, 20);
        WriteLn('First 20 digits: ', first20);
        WriteLn('Last 20 digits:  ', last20);
        WriteLn('Full number saved to memory');
    End;
    WriteLn;
End;

// ----------------------------------------------------------------------------
// Simple iterative Fibonacci (for comparison, limited to smaller n)
// ----------------------------------------------------------------------------
Procedure FibonacciIterative(n: Integer; Var result: mpz_t);
Var
    a, b, temp: mpz_t;
    i: Integer;
    abs_n: Integer;
    startTime, endTime: TDateTime;
Begin
    startTime := Now;

    mpz_init(a);
    mpz_init(b);
    mpz_init(temp);

    abs_n := Abs(n);

    If abs_n = 0 Then
    Begin
        mpz_set_si(result, 0);
    End
    Else
    Begin
        mpz_set_si(a, 0);
        mpz_set_si(b, 1);

        For i := 2 To abs_n Do
        Begin
            mpz_add(temp, a, b);
            mpz_set(a, b);
            mpz_set(b, temp);
        End;

        mpz_set(result, b);
    End;

    // Handle negative indices using: F(-n) = (-1)^(n+1) * F(n)
    If n < 0 Then
    Begin
        If (abs_n Mod 2) = 0 Then
            mpz_neg(result, result);
    End;

    mpz_clear(a);
    mpz_clear(b);
    mpz_clear(temp);

    endTime := Now;
    WriteLn('Iterative method time: ', MilliSecondsBetween(endTime, startTime), ' ms');
End;

// ----------------------------------------------------------------------------
// Parse a string to Int64, handling large numbers
// ----------------------------------------------------------------------------
Function ParseLargeInt(Str: String): Int64;
Var
    ErrorCode: Integer;
Begin
    Val(Str, Result, ErrorCode);
    If ErrorCode <> 0 Then
        Raise EConvertError.Create('Invalid integer: ' + Str);
End;

// ----------------------------------------------------------------------------
// Parse algorithm selection from string
// ----------------------------------------------------------------------------
Function ParseAlgorithm(Str: String): TAlgorithm;
Var
    lowerStr: String;
Begin
    lowerStr := LowerCase(Str);
    If lowerStr = 'matrix' Then
        Result := algMatrix
    Else If lowerStr = 'fastdoubling' Then
        Result := algFastDoubling
    Else If lowerStr = 'iterative' Then
        Result := algIterative
    Else If lowerStr = 'hybrid' Then
        Result := algHybrid
    Else
        Raise EConvertError.Create('Invalid algorithm: ' + Str + '. Use: matrix, fastdoubling, iterative, or hybrid');
End;

// ----------------------------------------------------------------------------
// Benchmark all methods for a given n
// ----------------------------------------------------------------------------
Procedure BenchmarkMethods(n: Int64);
Var
    result1, result2, result3, result4: mpz_t;
    s1, s2, s3, s4: AnsiString;
    startTime, endTime: TDateTime;
Begin
    WriteLn('========================================');
    WriteLn('Benchmarking F(', n, ')');
    WriteLn('========================================');

    mpz_init(result1);
    mpz_init(result2);
    mpz_init(result3);
    mpz_init(result4);

    // Benchmark Matrix method
    WriteLn;
    WriteLn('1. Matrix Exponentiation Method:');
    startTime := Now;
    FibonacciMatrix(n, result1);
    endTime := Now;
    s1 := mpz_get_str(Nil, 10, result1);
    WriteLn('Result length: ', Length(s1), ' digits');
    WriteLn('Total time: ', MilliSecondsBetween(endTime, startTime), ' ms');

    // Benchmark Fast Doubling method
    WriteLn;
    WriteLn('2. Fast Doubling Method:');
    startTime := Now;
    FibonacciFastDoubling(n, result2);
    endTime := Now;
    s2 := mpz_get_str(Nil, 10, result2);
    WriteLn('Result length: ', Length(s2), ' digits');
    WriteLn('Total time: ', MilliSecondsBetween(endTime, startTime), ' ms');

    // Benchmark Hybrid method
    WriteLn;
    WriteLn('3. Hybrid Doubling/Tripling Method:');
    startTime := Now;
    FibonacciHybrid(n, result3);
    endTime := Now;
    s3 := mpz_get_str(Nil, 10, result3);
    WriteLn('Result length: ', Length(s3), ' digits');
    WriteLn('Total time: ', MilliSecondsBetween(endTime, startTime), ' ms');

    // Verify results match
    If (mpz_cmp(result1, result2) = 0) and (mpz_cmp(result2, result3) = 0) Then
        WriteLn('✓ All results match!')
    Else
        WriteLn('✗ ERROR: Results do not match!');

    // Benchmark Iterative method if n is small enough
    If (n >= Low(Integer)) And (n <= High(Integer)) And (Abs(n) <= 100000) Then
    Begin
        WriteLn;
        WriteLn('4. Iterative Method (for comparison):');
        startTime := Now;
        FibonacciIterative(Integer(n), result4);
        endTime := Now;
        s4 := mpz_get_str(Nil, 10, result4);
        WriteLn('Result length: ', Length(s4), ' digits');
        WriteLn('Total time: ', MilliSecondsBetween(endTime, startTime), ' ms');

        If mpz_cmp(result1, result4) = 0 Then
            WriteLn('✓ Iterative result matches!')
        Else
            WriteLn('✗ ERROR: Iterative result does not match!');
    End
    Else
    Begin
        WriteLn;
        WriteLn('4. Iterative Method: Skipped (n too large)');
    End;

    WriteLn;

    mpz_clear(result1);
    mpz_clear(result2);
    mpz_clear(result3);
    mpz_clear(result4);
End;

// ----------------------------------------------------------------------------
// Main program with CLI argument handling
// ----------------------------------------------------------------------------
Var
    fib_gmp: mpz_t;
    n: Int64;
    i: Integer;
    algorithm: TAlgorithm;
    benchmarkMode: Boolean;

Begin
    // Initialize GMP result
    mpz_init(fib_gmp);

    // Check command line arguments
    If ParamCount >= 1 Then
    Begin
        // Parse first argument as 64-bit integer
        Try
            n := ParseLargeInt(ParamStr(1));

            // Check for algorithm selection (second argument)
            algorithm := algMatrix;  // Default
            benchmarkMode := False;

            For i := 2 To ParamCount Do
            Begin
                If LowerCase(ParamStr(i)) = '-benchmark' Then
                    benchmarkMode := True
                Else If Pos('-a=', LowerCase(ParamStr(i))) = 1 Then
                Begin
                    Try
                        algorithm := ParseAlgorithm(Trim(Copy(ParamStr(i), 4, MaxInt)));
                    Except
                        On E: EConvertError Do
                        Begin
                            WriteLn('Warning: ', E.Message);
                            WriteLn('Using default algorithm: matrix');
                        End;
                    End;
                End;
            End;

            If benchmarkMode Then
            Begin
                BenchmarkMethods(n);
            End
            Else
            Begin
                // Compute Fibonacci number with selected algorithm
                Case algorithm Of
                    algMatrix:
                    Begin
                        WriteLn('Computing F(', n, ') using matrix exponentiation...');
                        FibonacciMatrix(n, fib_gmp);
                    End;

                    algFastDoubling:
                    Begin
                        WriteLn('Computing F(', n, ') using fast doubling...');
                        FibonacciFastDoubling(n, fib_gmp);
                    End;

                    algIterative:
                    Begin
                        // For iterative method, we need to check if n fits in 32-bit
                        If (n < Low(Integer)) Or (n > High(Integer)) Then
                        Begin
                            WriteLn('Error: Iterative method only supports n between ',
                                    Low(Integer), ' and ', High(Integer));
                            WriteLn('Use matrix, fastdoubling, or hybrid method for larger values.');
                            Halt(1);
                        End;

                        WriteLn('Computing F(', n, ') using iterative method...');
                        FibonacciIterative(Integer(n), fib_gmp);
                    End;

                    algHybrid:
                    Begin
                        WriteLn('Computing F(', n, ') using hybrid (doubling/tripling) method...');
                        FibonacciHybrid(n, fib_gmp);
                    End;
                End;

                // Display result
                DisplayLargeFibonacci(n, fib_gmp);
            End;

        Except
            On E: EConvertError Do
                WriteLn('Error: ', E.Message);
        End;
    End
    Else
    Begin
        // No arguments: run demonstration
        WriteLn('No command line argument provided. Running demonstration...');
        WriteLn;
        WriteLn('========================================');
        WriteLn('Fibonacci Sequence Demonstration');
        WriteLn('========================================');
        WriteLn;

        // Test from -10 to 10
        WriteLn('Fibonacci numbers from -10 to 10 (using fast doubling):');
        For i := -10 To 10 Do
        Begin
            FibonacciFastDoubling(i, fib_gmp);
            WriteLn('F(', i:3, ') = ', mpz_get_str(Nil, 10, fib_gmp));
        End;
        WriteLn;

        // Benchmark different methods for moderate n
        WriteLn('Benchmark for n = 1,000,000:');
        BenchmarkMethods(1000000);

        // Benchmark for larger n
        WriteLn('Benchmark for n = 10,000,000:');
        BenchmarkMethods(10000000);

        // Comparison with iterative method for small n
        WriteLn('Comparison for small n (30):');
        Write('Matrix method: ');
        FibonacciMatrix(30, fib_gmp);
        WriteLn(mpz_get_str(Nil, 10, fib_gmp));

        Write('Hybrid method: ');
        FibonacciHybrid(30, fib_gmp);
        WriteLn(mpz_get_str(Nil, 10, fib_gmp));

        Write('Iterative method: ');
        FibonacciIterative(30, fib_gmp);
        WriteLn(mpz_get_str(Nil, 10, fib_gmp));
        WriteLn;

        WriteLn('Usage examples:');
        WriteLn('  ./FastFibonacciGMP 100                             # Compute F(100) with matrix (default)');
        WriteLn('  ./FastFibonacciGMP 100 -a=fastdoubling             # Compute with fast doubling');
        WriteLn('  ./FastFibonacciGMP 100 -a=hybrid                   # Compute with hybrid method (recommended)');
        WriteLn('  ./FastFibonacciGMP 100 -a=iterative                # Compute with iterative (small n only)');
        WriteLn('  ./FastFibonacciGMP 1000000 -benchmark              # Benchmark all methods');
        WriteLn('  ./FastFibonacciGMP -50 -a=fastdoubling             # Compute negative Fibonacci');
        WriteLn('  ./FastFibonacciGMP 18000000000 -a=hybrid # Compute large F(18^9) takes 11GB, 154701 ms to calculate/20 mins to write out');
        WriteLn;
        WriteLn('Available algorithms:');
        WriteLn('  matrix         : Matrix exponentiation (O(log n))');
        WriteLn('  fastdoubling   : Fast doubling iterative (O(log n), very fast)');
        WriteLn('  hybrid         : Hybrid doubling/tripling (O(log n), fastest for large n) ★ RECOMMENDED');
        WriteLn('  iterative      : Iterative method (O(n), for small n ≤ 100000 only)');
    End;

    // Cleanup
    mpz_clear(fib_gmp);

    WriteLn('Program completed.');
End.
(*)
*)  AI means corrected logic for AI itself: the first attempts were mistranslations
**) AI Finally understood what the algorithm is doing (after 12 sessions)
(*)
<PRE>