Skip to content

Instantly share code, notes, and snippets.

@Thraetaona
Last active November 14, 2024 17:23
Show Gist options
  • Save Thraetaona/1597b89e28980b26f156ce50ece40bb3 to your computer and use it in GitHub Desktop.
Save Thraetaona/1597b89e28980b26f156ce50ece40bb3 to your computer and use it in GitHub Desktop.
Newton's Divided Differences table generator implemented in the Ada language.
-- ---------------------------------------------------------------------
-- 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
-- ---------------------------------------------------------------------
@Thraetaona
Copy link
Author

Thraetaona commented Oct 31, 2024

-- ---------------------------------------------------------------------
-- 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
-- ---------------------------------------------------------------------

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment