Newton's Divided Differences table generator implemented in the Ada language.

divided_differences.adb

-- ---------------------------------------------------------------------
-- SPDX-License-Identifier: CC0-1.0
--
-- Authored by Fereydoun Memarzanjany
--
-- To the maximum extent possible under law, the Author waives all
-- copyright and related or neighboring rights to this code.
--
-- You should have received a copy of the CC0 legalcode along with this
-- work; if not, see <http://creativecommons.org/publicdomain/zero/1.0/>
-- ---------------------------------------------------------------------


with Ada.Text_IO,
     Ada.Integer_Text_IO,
     Ada.Float_Text_IO;
    use Ada.Text_IO,
        Ada.Integer_Text_IO,
        Ada.Float_Text_IO;

with Ada.Command_Line;
    use Ada.Command_Line;

procedure Divided_Differences_Table is
    n : Integer;
begin
    Read_And_Validate_Input : declare
    begin
        <<Get_Input>>
        Put("Enter the number of data points (n): ");
        Get(n);
        Skip_Line;

        <<Validate_Input>>
        if n < 1 then
            Put_Line("Invalid number of data points.");
            Set_Exit_Status(Failure);
            return;
        end if;
    end Read_And_Validate_Input;

    -- -----------------------------------------------------------------
    -- Initialize the table and compute divided differences
    -- -----------------------------------------------------------------

    Setup_Table : declare
        x    : array (1 .. n) of Float;
        f    : array (1 .. n) of Float;
        -- Divided differences table
        diff : array (1 .. n, 0 .. n - 1) of Float;
        i, j : Positive;
    begin

        -- -------------------------------------------------------------
        -- Obtain values from the user
        -- -------------------------------------------------------------
        <<Read_Values>>

        Put_Line("Enter x and f(x) values:");
        for i in 1 .. n loop
            Put("x("); Put(i, Width => 0); Put(") = ");
            Get(x(i));
            Skip_Line;
            Put("f("); Put(i, Width => 0); Put(") = ");
            Get(f(i));
            Skip_Line;
        end loop;

        <<Initialize_Table>>
        -- Initialize the 1st dimension (i.e., rows) with f(x) values
        for i in diff'Range(1) loop
            diff(i, 0) := f(i);
        end loop;

        -- -------------------------------------------------------------
        -- Compute divided differences
        -- -------------------------------------------------------------
        <<Compute_Differences>>

        -- Iterate over columns (orders of divided differences)
        for j in diff'First(2) + 1 .. diff'Last(2) loop
            -- Iterate over rows for the current order/column
            for i in diff'First(1) .. diff'Last(1) - j loop
                diff(i, j) :=
                    (diff(i + 1, j - 1) - diff(i, j - 1)) /
                    (x(i + j) - x(i));
            end loop;
        end loop;

        -- -------------------------------------------------------------
        -- Display the table
        -- -------------------------------------------------------------
        <<Display_Output>>

        Print_Table_Header : declare
        begin
            -- Print the column headers of the divided differences table
            -- (This first prints i, x(i), f(xi) and then prints
            -- "Order <N>" for each higher order difference.)
            New_Line;
            Put_Line("Divided Differences Table:");
            Put("i   x(i)       f(xi)      ");
            -- Print "Order <N>" per each column now.
            for j in diff'First(2) + 1 .. diff'Last(2) loop
                Put("Order "); Put(j, Width => 0); Put("       ");
            end loop;
            New_Line;
        end Print_Table_Header;

        Print_Table_Data : declare
            -- Local procedure for consistent float formatting
            procedure Put_Float(Item : Float) is
            begin
                Put(Item, Fore => 8, Aft => 5, Exp => 0);
            end Put_Float;
        begin
            -- Print each row of the divided differences table.
            -- (This first prints row index i, x value, f(x) value, and
            -- then it prints the divided differences for each order.)
                for i in diff'Range(1) loop
                    Put(i, Width => 2); Put("  ");
                    Put_Float(x(i)); Put("  ");
                    Put_Float(diff(i, 0));

                    for j in diff'First(2) + 1 ..
                        diff'Last(2) - (i - diff'First(1))
                    loop
                        Put("  ");
                        Put_Float(diff(i, j));
                    end loop;
                    New_Line;
                end loop;
        end Print_Table_Data;

    Compute_Lagrange : declare
        num_points : Integer;
        x_eval     : array (1 .. 100) of Float;  -- Array for evaluation 
        result     : Float;                      -- Result of L(x)
        term       : Float;                      -- Temporary storage
    begin
        -- Get number of evaluation points
        Put("Enter number of x values to evaluate: ");
        Get(num_points);
        Skip_Line;

        -- Get all evaluation points
        for i in 1 .. num_points loop
            Put("Enter x("); Put(i, Width => 0);
            Put(") to evaluate L(x) at: ");
            Get(x_eval(i));
            Skip_Line;
        end loop;

        -- Calculate and display results for each point
        New_Line;
        Put_Line("Results:");
        for i in 1 .. num_points loop
            -- Start with the constant term
            result := diff(1, 0);

            -- Add each term in the Newton form
            for j in 1 .. n-1 loop
                -- Compute (x - x0)(x - x1)...(x - x_(j-1))
                term := diff(1, j);
                for k in 1 .. j loop
                    term := term * (x_eval(i) - x(k));
                end loop;
                result := result + term;
            end loop;

            -- Display result
            Put("L(");
            Put(x_eval(i), Fore => 0, Aft => 5, Exp => 0);
            Put(") = ");
            Put(result, Fore => 0, Aft => 5, Exp => 0);
            New_Line;
        end loop;
    end Compute_Lagrange;

    end Setup_Table;
end Divided_Differences_Table;


-- ---------------------------------------------------------------------
-- END OF FILE: main.adb
-- ---------------------------------------------------------------------

-- ---------------------------------------------------------------------
-- Sample Output:
-- ---------------------------------------------------------------------
-- Enter the number of data points (n): 4
-- Enter x and f(x) values:
-- x(1) = 0
-- f(1) = 1
-- x(2) = 1
-- f(2) = 2.718282
-- x(3) = 1.5
-- f(3) = 4.481689
-- x(4) = 2
-- f(4) = 7.389056
-- 
-- Divided Differences Table:
-- i  x(i)     f(xi)    Order 1  Order 2  Order 3
--  1  0.00000  1.00000  1.71828  1.20569  0.54112
--  2  1.00000  2.71828  3.52681  2.28792
--  3  1.50000  4.48169  5.81473
--  4  2.00000  7.38906
--
-- Enter number of x values to evaluate: 8
-- Enter x(1) to evaluate L(x) at: -1
-- Enter x(2) to evaluate L(x) at: 0
-- Enter x(3) to evaluate L(x) at: 0.5
-- Enter x(4) to evaluate L(x) at: 0.75
-- Enter x(5) to evaluate L(x) at: 1
-- Enter x(6) to evaluate L(x) at: 1.5
-- Enter x(8) to evaluate L(x) at: 15
--
-- Results:
-- L(-1.00000) = -1.01249
-- L(0.00000) = 1.00000
-- L(0.50000) = 1.69300
-- L(0.75000) = 2.13874
-- L(1.00000) = 2.71828
-- L(1.50000) = 4.48169
-- L(2.00000) = 7.38906
-- ---------------------------------------------------------------------

-- ---------------------------------------------------------------------
-- Sample Output:
-- ---------------------------------------------------------------------
-- Enter the number of data points (n): 4
-- Enter x and f(x) values:
-- x(1) = 0
-- f(1) = 1
-- x(2) = 1
-- f(2) = 2.718282
-- x(3) = 1.5
-- f(3) = 4.481689
-- x(4) = 2
-- f(4) = 7.389056
-- 
-- Divided Differences Table:
-- i  x(i)     f(xi)    Order 1  Order 2  Order 3
--  1  0.00000  1.00000  1.71828  1.20569  0.54112
--  2  1.00000  2.71828  3.52681  2.28792
--  3  1.50000  4.48169  5.81473
--  4  2.00000  7.38906
--
-- Enter number of x values to evaluate: 8
-- Enter x(1) to evaluate L(x) at: -1
-- Enter x(2) to evaluate L(x) at: 0
-- Enter x(3) to evaluate L(x) at: 0.5
-- Enter x(4) to evaluate L(x) at: 0.75
-- Enter x(5) to evaluate L(x) at: 1
-- Enter x(6) to evaluate L(x) at: 1.5
-- Enter x(8) to evaluate L(x) at: 15
--
-- Results:
-- L(-1.00000) = -1.01249
-- L(0.00000) = 1.00000
-- L(0.50000) = 1.69300
-- L(0.75000) = 2.13874
-- L(1.00000) = 2.71828
-- L(1.50000) = 4.48169
-- L(2.00000) = 7.38906
-- ---------------------------------------------------------------------