1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147 | with System;
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Text_IO.C_Streams;
with Ada.Long_Float_Text_IO; use Ada.Long_Float_Text_IO;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
with Interfaces.C; use Interfaces.C;
with Interfaces.C.Strings; use Interfaces.C.Strings;
with Ada.Strings.Unbounded; use Ada.Strings.Unbounded;
with Ada.Command_Line; use Ada.Command_Line;
with Ada.Numerics.Elementary_Functions ; use Ada.Numerics.Elementary_Functions ;
with GNAT.Source_Info; use GNAT.Source_Info;
with gsl;
with gsl.math ;
with gsl.integration ;
with gsl.sf.gamma ;
with fns ;
procedure integration is
use gsl.double_text_Io;
use gsl.double_elementary_functions ;
Verbose : boolean := True ;
Status : Int ;
-- 17.16.1 Adaptive integration example
procedure Test1 is
myname : String := gnat.Source_Info.enclosing_entity ;
logfile : File_Type ;
N : constant := 1000 ;
ws : access gsl.integration.gsl_integration_workspace
:= gsl.integration.workspace_alloc(N) ;
result, error : aliased double ;
expected : constant double := -4.0;
alpha : aliased double := 1.0;
fn : aliased gsl.math.gsl_function ;
begin
if Verbose
then
Put_Line(myname);
end if ;
fn.c_function := fns.LogSqrt'Access ;
fn.params := alpha'Address ;
Status := gsl.integration.qags( fn'Access , 0.0 , 1.0 , 0.0 , 1.0e-7 , 1000, ws, result'access , error'access);
Put("Result "); Put(result, aft => 6 ) ; Put(" ");
Put("Rel Error "); Put(error) ; Put(" ");
Put("Intervals "); Put(Integer(ws.size));
New_Line;
-- Create(logfile,Out_File,myname & ".csv");
-- Set_Output(logfile);
-- Set_Output(Standard_Output);
-- Close(logfile);
gsl.integration.workspace_free(ws);
end Test1 ;
-- 17.16.2 Fixed-point quadrature example
procedure Test2 is
myname : String := gnat.Source_Info.enclosing_entity ;
logfile : File_Type ;
N : constant := 6 ;
M : aliased int := 10 ;
ws : access gsl.integration.gsl_integration_fixed_workspace ;
fn : aliased gsl.math.gsl_function ;
result : aliased double ;
begin
if Verbose
then
Put_Line(myname);
end if ;
ws := gsl.integration.fixed_alloc(gsl.integration.gsl_integration_fixed_hermite,
N , 0.0 , 1.0 , 0.0 , 0.0 );
fn.c_function := fns.Power'access ;
fn.params := M'Address ;
Status := gsl.integration.fixed(fn'access, result'access, ws ) ;
Put("Result "); Put(result); New_Line ;
Put("Expected ");
if m mod 2 = 0
then
Put(gsl.math.M_SQRTPI + gsl.sf.gamma.gamma(0.5*(1.0 + double(m))));
else
Put(double(gsl.math.M_SQRTPI));
end if ;
New_Line;
Create(logfile,Out_File,myname & ".csv");
Set_Output(logfile);
Set_Output(Standard_Output);
Close(logfile);
gsl.integration.fixed_free(ws);
end Test2 ;
procedure Test3 is
myname : String := gnat.Source_Info.enclosing_entity ;
logfile : File_Type ;
N : constant := 6 ;
ws : access gsl.integration.gsl_integration_fixed_workspace ;
fn : aliased gsl.math.gsl_function ;
result : aliased double ;
x : double ;
xdelta : constant double := 2.0 * Ada.Numerics.Pi/256.0 ;
i : Integer := 1 ;
begin
if Verbose
then
Put_Line(myname);
end if ;
fn.c_function := fns.Sin'access ;
fn.params := System.Null_Address ;
Create(logfile,Out_File,myname & ".csv");
Set_Output(logfile);
x := xdelta;
i := 0 ;
while i < 256
loop
i := i + 1 ;
ws := gsl.integration.fixed_alloc(gsl.integration.gsl_integration_fixed_legendre,
N , 0.0 , x , 0.0 , 0.0 );
Put(i) ; Put(" ; ");
Put(x) ; Put(" ; ");
Put(fns.Sin(x,System.Null_Address)); Put(" ; ");
Status := gsl.integration.fixed(fn'access , result'access , ws);
Put(result) ; Put(" ; ");
New_Line ;
x := x + xdelta ;
gsl.integration.fixed_free(ws);
end loop ;
Set_Output(Standard_Output);
Close(logfile);
end Test3 ;
begin
if Verbose
then
Put_line(Gnat.Source_Info.enclosing_entity);
end if;
gsl.double_text_Io.Default_exp := 0;
Test1 ;
Test2 ;
Test3 ;
end integration;
|