gslshell/integ.lua at master · LuaDist/gslshell · GitHub

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
-- integr module follows

localfunctionget_spec(spec, key, mtype)
localdefaults= {eps_abs=1e-6, eps_rel=1e-6, limit=512,
rule='SINGULAR'}
localv=spec[key]
ifvandtype(v) ~=mtypethen
error(string.format('parameter %s should be of type %s', key, mtype))
end
returnvandvordefaults[key]
end

localfunctiontype_check(v, mtype, name)
ifnotvthenerror('parameter ' ..name..' required') end
iftype(v) ~=mtypethen
error('parameter ' ..name..' should be of type ' ..mtype)
end
end

localinteg_raw=integ_module_init ()

functioninteg(spec)
locallimit_min=512

localf=spec.f
type_check (spec.f, 'function', 'f')
ifnotspec.pointsor#spec.points<2then
error('integration points not given')
end
localpts=spec.points
localinttype

localcspec= {}
localinteg_spec_add=
function(key, mtype)
localvalue=get_spec(spec, key, mtype)
cspec[key] =value
end

integ_spec_add ('eps_abs', 'number')
integ_spec_add ('eps_rel', 'number')
integ_spec_add ('limit', 'number')
ifcspec.limit<limit_minthencspec.limit=limit_minend

cspec.a=pts[1]
cspec.b=pts[#pts]

-- infinite intervals
ifpts[1] =='-inf' orpts[#pts] =='+inf' then
ifpts[1] =='-inf' andpts[#pts] =='+inf' then
inttype='agi'
elseifpts[1] =='-inf' then
inttype='agil'
else
inttype='agiu'
end
ifspec.weightthen
localw=spec.weight
ifw.type=='sin' orw.type=='cos' then
inttype=string.gsub(inttype, '^agi', 'awf')
if (inttype=='awf') then
cspec.eps_abs=cspec.eps_abs/2
cspec.a, cspec.b=0, 0
locallint, lerr=integ_raw(f, 'awfl', cspec, w)
localuint, uerr=integ_raw(f, 'awfu', cspec, w)
returnlint+uint, lerr+uerr
end
returninteg_raw(f, inttype, cspec, w)
end
error('cannot calculate indefinite integral with this weight')
end
returninteg_raw(f, inttype, cspec)
end

fori,vinipairs(pts) do
iftype(v) ~='number' thenerror('invalid points') end
end

-- weighted integrals
ifspec.weightthen
localw=spec.weight
if#pts>2then
error('weighted integrals does not accept singular points')
end
ifw.type=='sin' orw.type=='cos' then
type_check (w.omega, 'number', 'omega')
inttype='awo'
elseifw.type=='cauchy' then
type_check (w.singularity, 'number', 'singularity')
inttype='awc'
elseifw.type=='alglog' then
inttype='aws'
else
error('unknown wight type')
end
returninteg_raw(f, inttype, cspec, w)
end

ifnotspec.adaptivethen
if#pts>2then
error('non-adaptive integrals does not accept singular points')
end
inttype='ng'
end

if#pts>2then
integ_spec_add ('points', 'table')
inttype='agp'
end

localrulename=get_spec(spec, 'rule', 'string')
ifrulename=='SINGULAR' then
inttype='ags'
else
cspec.rule=tonumber(string.match(rulename, 'GAUSS(%d+)'))
ifnotcspec.rulethenerror('invalid integration rule') end
inttype='ag'
end

returninteg_raw(f, inttype, cspec)
end
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
-- integr module follows

localfunctionget_spec(spec, key, mtype)
localdefaults= {eps_abs=1e-6, eps_rel=1e-6, limit=512,
rule='SINGULAR'}
localv=spec[key]
ifvandtype(v) ~=mtypethen
error(string.format('parameter %s should be of type %s', key, mtype))
end
returnvandvordefaults[key]
end

localfunctiontype_check(v, mtype, name)
ifnotvthenerror('parameter ' ..name..' required') end
iftype(v) ~=mtypethen
error('parameter ' ..name..' should be of type ' ..mtype)
end
end

localinteg_raw=integ_module_init ()

functioninteg(spec)
locallimit_min=512

localf=spec.f
type_check (spec.f, 'function', 'f')
ifnotspec.pointsor#spec.points<2then
error('integration points not given')
end
localpts=spec.points
localinttype

localcspec= {}
localinteg_spec_add=
function(key, mtype)
localvalue=get_spec(spec, key, mtype)
cspec[key] =value
end

integ_spec_add ('eps_abs', 'number')
integ_spec_add ('eps_rel', 'number')
integ_spec_add ('limit', 'number')
ifcspec.limit<limit_minthencspec.limit=limit_minend

cspec.a=pts[1]
cspec.b=pts[#pts]

-- infinite intervals
ifpts[1] =='-inf' orpts[#pts] =='+inf' then
ifpts[1] =='-inf' andpts[#pts] =='+inf' then
inttype='agi'
elseifpts[1] =='-inf' then
inttype='agil'
else
inttype='agiu'
end
ifspec.weightthen
localw=spec.weight
ifw.type=='sin' orw.type=='cos' then
inttype=string.gsub(inttype, '^agi', 'awf')
if (inttype=='awf') then
cspec.eps_abs=cspec.eps_abs/2
cspec.a, cspec.b=0, 0
locallint, lerr=integ_raw(f, 'awfl', cspec, w)
localuint, uerr=integ_raw(f, 'awfu', cspec, w)
returnlint+uint, lerr+uerr
end
returninteg_raw(f, inttype, cspec, w)
end
error('cannot calculate indefinite integral with this weight')
end
returninteg_raw(f, inttype, cspec)
end

fori,vinipairs(pts) do
iftype(v) ~='number' thenerror('invalid points') end
end

-- weighted integrals
ifspec.weightthen
localw=spec.weight
if#pts>2then
error('weighted integrals does not accept singular points')
end
ifw.type=='sin' orw.type=='cos' then
type_check (w.omega, 'number', 'omega')
inttype='awo'
elseifw.type=='cauchy' then
type_check (w.singularity, 'number', 'singularity')
inttype='awc'
elseifw.type=='alglog' then
inttype='aws'
else
error('unknown wight type')
end
returninteg_raw(f, inttype, cspec, w)
end

ifnotspec.adaptivethen
if#pts>2then
error('non-adaptive integrals does not accept singular points')
end
inttype='ng'
end

if#pts>2then
integ_spec_add ('points', 'table')
inttype='agp'
end

localrulename=get_spec(spec, 'rule', 'string')
ifrulename=='SINGULAR' then
inttype='ags'
else
cspec.rule=tonumber(string.match(rulename, 'GAUSS(%d+)'))
ifnotcspec.rulethenerror('invalid integration rule') end
inttype='ag'
end

returninteg_raw(f, inttype, cspec)
end
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
-- integr module follows

localfunctionget_spec(spec, key, mtype)
localdefaults= {eps_abs=1e-6, eps_rel=1e-6, limit=512,
rule='SINGULAR'}
localv=spec[key]
ifvandtype(v) ~=mtypethen
error(string.format('parameter %s should be of type %s', key, mtype))
end
returnvandvordefaults[key]
end

localfunctiontype_check(v, mtype, name)
ifnotvthenerror('parameter ' ..name..' required') end
iftype(v) ~=mtypethen
error('parameter ' ..name..' should be of type ' ..mtype)
end
end

localinteg_raw=integ_module_init ()

functioninteg(spec)
locallimit_min=512

localf=spec.f
type_check (spec.f, 'function', 'f')
ifnotspec.pointsor#spec.points<2then
error('integration points not given')
end
localpts=spec.points
localinttype

localcspec= {}
localinteg_spec_add=
function(key, mtype)
localvalue=get_spec(spec, key, mtype)
cspec[key] =value
end

integ_spec_add ('eps_abs', 'number')
integ_spec_add ('eps_rel', 'number')
integ_spec_add ('limit', 'number')
ifcspec.limit<limit_minthencspec.limit=limit_minend

cspec.a=pts[1]
cspec.b=pts[#pts]

-- infinite intervals
ifpts[1] =='-inf' orpts[#pts] =='+inf' then
ifpts[1] =='-inf' andpts[#pts] =='+inf' then
inttype='agi'
elseifpts[1] =='-inf' then
inttype='agil'
else
inttype='agiu'
end
ifspec.weightthen
localw=spec.weight
ifw.type=='sin' orw.type=='cos' then
inttype=string.gsub(inttype, '^agi', 'awf')
if (inttype=='awf') then
cspec.eps_abs=cspec.eps_abs/2
cspec.a, cspec.b=0, 0
locallint, lerr=integ_raw(f, 'awfl', cspec, w)
localuint, uerr=integ_raw(f, 'awfu', cspec, w)
returnlint+uint, lerr+uerr
end
returninteg_raw(f, inttype, cspec, w)
end
error('cannot calculate indefinite integral with this weight')
end
returninteg_raw(f, inttype, cspec)
end

fori,vinipairs(pts) do
iftype(v) ~='number' thenerror('invalid points') end
end

-- weighted integrals
ifspec.weightthen
localw=spec.weight
if#pts>2then
error('weighted integrals does not accept singular points')
end
ifw.type=='sin' orw.type=='cos' then
type_check (w.omega, 'number', 'omega')
inttype='awo'
elseifw.type=='cauchy' then
type_check (w.singularity, 'number', 'singularity')
inttype='awc'
elseifw.type=='alglog' then
inttype='aws'
else
error('unknown wight type')
end
returninteg_raw(f, inttype, cspec, w)
end

ifnotspec.adaptivethen
if#pts>2then
error('non-adaptive integrals does not accept singular points')
end
inttype='ng'
end

if#pts>2then
integ_spec_add ('points', 'table')
inttype='agp'
end

localrulename=get_spec(spec, 'rule', 'string')
ifrulename=='SINGULAR' then
inttype='ags'
else
cspec.rule=tonumber(string.match(rulename, 'GAUSS(%d+)'))
ifnotcspec.rulethenerror('invalid integration rule') end
inttype='ag'
end

returninteg_raw(f, inttype, cspec)
end
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
-- integr module follows

localfunctionget_spec(spec, key, mtype)
localdefaults= {eps_abs=1e-6, eps_rel=1e-6, limit=512,
rule='SINGULAR'}
localv=spec[key]
ifvandtype(v) ~=mtypethen
error(string.format('parameter %s should be of type %s', key, mtype))
end
returnvandvordefaults[key]
end

localfunctiontype_check(v, mtype, name)
ifnotvthenerror('parameter ' ..name..' required') end
iftype(v) ~=mtypethen
error('parameter ' ..name..' should be of type ' ..mtype)
end
end

localinteg_raw=integ_module_init ()

functioninteg(spec)
locallimit_min=512

localf=spec.f
type_check (spec.f, 'function', 'f')
ifnotspec.pointsor#spec.points<2then
error('integration points not given')
end
localpts=spec.points
localinttype

localcspec= {}
localinteg_spec_add=
function(key, mtype)
localvalue=get_spec(spec, key, mtype)
cspec[key] =value
end

integ_spec_add ('eps_abs', 'number')
integ_spec_add ('eps_rel', 'number')
integ_spec_add ('limit', 'number')
ifcspec.limit<limit_minthencspec.limit=limit_minend

cspec.a=pts[1]
cspec.b=pts[#pts]

-- infinite intervals
ifpts[1] =='-inf' orpts[#pts] =='+inf' then
ifpts[1] =='-inf' andpts[#pts] =='+inf' then
inttype='agi'
elseifpts[1] =='-inf' then
inttype='agil'
else
inttype='agiu'
end
ifspec.weightthen
localw=spec.weight
ifw.type=='sin' orw.type=='cos' then
inttype=string.gsub(inttype, '^agi', 'awf')
if (inttype=='awf') then
cspec.eps_abs=cspec.eps_abs/2
cspec.a, cspec.b=0, 0
locallint, lerr=integ_raw(f, 'awfl', cspec, w)
localuint, uerr=integ_raw(f, 'awfu', cspec, w)
returnlint+uint, lerr+uerr
end
returninteg_raw(f, inttype, cspec, w)
end
error('cannot calculate indefinite integral with this weight')
end
returninteg_raw(f, inttype, cspec)
end

fori,vinipairs(pts) do
iftype(v) ~='number' thenerror('invalid points') end
end

-- weighted integrals
ifspec.weightthen
localw=spec.weight
if#pts>2then
error('weighted integrals does not accept singular points')
end
ifw.type=='sin' orw.type=='cos' then
type_check (w.omega, 'number', 'omega')
inttype='awo'
elseifw.type=='cauchy' then
type_check (w.singularity, 'number', 'singularity')
inttype='awc'
elseifw.type=='alglog' then
inttype='aws'
else
error('unknown wight type')
end
returninteg_raw(f, inttype, cspec, w)
end

ifnotspec.adaptivethen
if#pts>2then
error('non-adaptive integrals does not accept singular points')
end
inttype='ng'
end

if#pts>2then
integ_spec_add ('points', 'table')
inttype='agp'
end

localrulename=get_spec(spec, 'rule', 'string')
ifrulename=='SINGULAR' then
inttype='ags'
else
cspec.rule=tonumber(string.match(rulename, 'GAUSS(%d+)'))
ifnotcspec.rulethenerror('invalid integration rule') end
inttype='ag'
end

returninteg_raw(f, inttype, cspec)
end
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
-- integr module follows

localfunctionget_spec(spec, key, mtype)
localdefaults= {eps_abs=1e-6, eps_rel=1e-6, limit=512,
rule='SINGULAR'}
localv=spec[key]
ifvandtype(v) ~=mtypethen
error(string.format('parameter %s should be of type %s', key, mtype))
end
returnvandvordefaults[key]
end

localfunctiontype_check(v, mtype, name)
ifnotvthenerror('parameter ' ..name..' required') end
iftype(v) ~=mtypethen
error('parameter ' ..name..' should be of type ' ..mtype)
end
end

localinteg_raw=integ_module_init ()

functioninteg(spec)
locallimit_min=512

localf=spec.f
type_check (spec.f, 'function', 'f')
ifnotspec.pointsor#spec.points<2then
error('integration points not given')
end
localpts=spec.points
localinttype

localcspec= {}
localinteg_spec_add=
function(key, mtype)
localvalue=get_spec(spec, key, mtype)
cspec[key] =value
end

integ_spec_add ('eps_abs', 'number')
integ_spec_add ('eps_rel', 'number')
integ_spec_add ('limit', 'number')
ifcspec.limit<limit_minthencspec.limit=limit_minend

cspec.a=pts[1]
cspec.b=pts[#pts]

-- infinite intervals
ifpts[1] =='-inf' orpts[#pts] =='+inf' then
ifpts[1] =='-inf' andpts[#pts] =='+inf' then
inttype='agi'
elseifpts[1] =='-inf' then
inttype='agil'
else
inttype='agiu'
end
ifspec.weightthen
localw=spec.weight
ifw.type=='sin' orw.type=='cos' then
inttype=string.gsub(inttype, '^agi', 'awf')
if (inttype=='awf') then
cspec.eps_abs=cspec.eps_abs/2
cspec.a, cspec.b=0, 0
locallint, lerr=integ_raw(f, 'awfl', cspec, w)
localuint, uerr=integ_raw(f, 'awfu', cspec, w)
returnlint+uint, lerr+uerr
end
returninteg_raw(f, inttype, cspec, w)
end
error('cannot calculate indefinite integral with this weight')
end
returninteg_raw(f, inttype, cspec)
end

fori,vinipairs(pts) do
iftype(v) ~='number' thenerror('invalid points') end
end

-- weighted integrals
ifspec.weightthen
localw=spec.weight
if#pts>2then
error('weighted integrals does not accept singular points')
end
ifw.type=='sin' orw.type=='cos' then
type_check (w.omega, 'number', 'omega')
inttype='awo'
elseifw.type=='cauchy' then
type_check (w.singularity, 'number', 'singularity')
inttype='awc'
elseifw.type=='alglog' then
inttype='aws'
else
error('unknown wight type')
end
returninteg_raw(f, inttype, cspec, w)
end

ifnotspec.adaptivethen
if#pts>2then
error('non-adaptive integrals does not accept singular points')
end
inttype='ng'
end

if#pts>2then
integ_spec_add ('points', 'table')
inttype='agp'
end

localrulename=get_spec(spec, 'rule', 'string')
ifrulename=='SINGULAR' then
inttype='ags'
else
cspec.rule=tonumber(string.match(rulename, 'GAUSS(%d+)'))
ifnotcspec.rulethenerror('invalid integration rule') end
inttype='ag'
end

returninteg_raw(f, inttype, cspec)
end
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
-- integr module follows

localfunctionget_spec(spec, key, mtype)
localdefaults= {eps_abs=1e-6, eps_rel=1e-6, limit=512,
rule='SINGULAR'}
localv=spec[key]
ifvandtype(v) ~=mtypethen
error(string.format('parameter %s should be of type %s', key, mtype))
end
returnvandvordefaults[key]
end

localfunctiontype_check(v, mtype, name)
ifnotvthenerror('parameter ' ..name..' required') end
iftype(v) ~=mtypethen
error('parameter ' ..name..' should be of type ' ..mtype)
end
end

localinteg_raw=integ_module_init ()

functioninteg(spec)
locallimit_min=512

localf=spec.f
type_check (spec.f, 'function', 'f')
ifnotspec.pointsor#spec.points<2then
error('integration points not given')
end
localpts=spec.points
localinttype

localcspec= {}
localinteg_spec_add=
function(key, mtype)
localvalue=get_spec(spec, key, mtype)
cspec[key] =value
end

integ_spec_add ('eps_abs', 'number')
integ_spec_add ('eps_rel', 'number')
integ_spec_add ('limit', 'number')
ifcspec.limit<limit_minthencspec.limit=limit_minend

cspec.a=pts[1]
cspec.b=pts[#pts]

-- infinite intervals
ifpts[1] =='-inf' orpts[#pts] =='+inf' then
ifpts[1] =='-inf' andpts[#pts] =='+inf' then
inttype='agi'
elseifpts[1] =='-inf' then
inttype='agil'
else
inttype='agiu'
end
ifspec.weightthen
localw=spec.weight
ifw.type=='sin' orw.type=='cos' then
inttype=string.gsub(inttype, '^agi', 'awf')
if (inttype=='awf') then
cspec.eps_abs=cspec.eps_abs/2
cspec.a, cspec.b=0, 0
locallint, lerr=integ_raw(f, 'awfl', cspec, w)
localuint, uerr=integ_raw(f, 'awfu', cspec, w)
returnlint+uint, lerr+uerr
end
returninteg_raw(f, inttype, cspec, w)
end
error('cannot calculate indefinite integral with this weight')
end
returninteg_raw(f, inttype, cspec)
end

fori,vinipairs(pts) do
iftype(v) ~='number' thenerror('invalid points') end
end

-- weighted integrals
ifspec.weightthen
localw=spec.weight
if#pts>2then
error('weighted integrals does not accept singular points')
end
ifw.type=='sin' orw.type=='cos' then
type_check (w.omega, 'number', 'omega')
inttype='awo'
elseifw.type=='cauchy' then
type_check (w.singularity, 'number', 'singularity')
inttype='awc'
elseifw.type=='alglog' then
inttype='aws'
else
error('unknown wight type')
end
returninteg_raw(f, inttype, cspec, w)
end

ifnotspec.adaptivethen
if#pts>2then
error('non-adaptive integrals does not accept singular points')
end
inttype='ng'
end

if#pts>2then
integ_spec_add ('points', 'table')
inttype='agp'
end

localrulename=get_spec(spec, 'rule', 'string')
ifrulename=='SINGULAR' then
inttype='ags'
else
cspec.rule=tonumber(string.match(rulename, 'GAUSS(%d+)'))
ifnotcspec.rulethenerror('invalid integration rule') end
inttype='ag'
end

returninteg_raw(f, inttype, cspec)
end
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
-- integr module follows

localfunctionget_spec(spec, key, mtype)
localdefaults= {eps_abs=1e-6, eps_rel=1e-6, limit=512,
rule='SINGULAR'}
localv=spec[key]
ifvandtype(v) ~=mtypethen
error(string.format('parameter %s should be of type %s', key, mtype))
end
returnvandvordefaults[key]
end

localfunctiontype_check(v, mtype, name)
ifnotvthenerror('parameter ' ..name..' required') end
iftype(v) ~=mtypethen
error('parameter ' ..name..' should be of type ' ..mtype)
end
end

localinteg_raw=integ_module_init ()

functioninteg(spec)
locallimit_min=512

localf=spec.f
type_check (spec.f, 'function', 'f')
ifnotspec.pointsor#spec.points<2then
error('integration points not given')
end
localpts=spec.points
localinttype

localcspec= {}
localinteg_spec_add=
function(key, mtype)
localvalue=get_spec(spec, key, mtype)
cspec[key] =value
end

integ_spec_add ('eps_abs', 'number')
integ_spec_add ('eps_rel', 'number')
integ_spec_add ('limit', 'number')
ifcspec.limit<limit_minthencspec.limit=limit_minend

cspec.a=pts[1]
cspec.b=pts[#pts]

-- infinite intervals
ifpts[1] =='-inf' orpts[#pts] =='+inf' then
ifpts[1] =='-inf' andpts[#pts] =='+inf' then
inttype='agi'
elseifpts[1] =='-inf' then
inttype='agil'
else
inttype='agiu'
end
ifspec.weightthen
localw=spec.weight
ifw.type=='sin' orw.type=='cos' then
inttype=string.gsub(inttype, '^agi', 'awf')
if (inttype=='awf') then
cspec.eps_abs=cspec.eps_abs/2
cspec.a, cspec.b=0, 0
locallint, lerr=integ_raw(f, 'awfl', cspec, w)
localuint, uerr=integ_raw(f, 'awfu', cspec, w)
returnlint+uint, lerr+uerr
end
returninteg_raw(f, inttype, cspec, w)
end
error('cannot calculate indefinite integral with this weight')
end
returninteg_raw(f, inttype, cspec)
end

fori,vinipairs(pts) do
iftype(v) ~='number' thenerror('invalid points') end
end

-- weighted integrals
ifspec.weightthen
localw=spec.weight
if#pts>2then
error('weighted integrals does not accept singular points')
end
ifw.type=='sin' orw.type=='cos' then
type_check (w.omega, 'number', 'omega')
inttype='awo'
elseifw.type=='cauchy' then
type_check (w.singularity, 'number', 'singularity')
inttype='awc'
elseifw.type=='alglog' then
inttype='aws'
else
error('unknown wight type')
end
returninteg_raw(f, inttype, cspec, w)
end

ifnotspec.adaptivethen
if#pts>2then
error('non-adaptive integrals does not accept singular points')
end
inttype='ng'
end

if#pts>2then
integ_spec_add ('points', 'table')
inttype='agp'
end

localrulename=get_spec(spec, 'rule', 'string')
ifrulename=='SINGULAR' then
inttype='ags'
else
cspec.rule=tonumber(string.match(rulename, 'GAUSS(%d+)'))
ifnotcspec.rulethenerror('invalid integration rule') end
inttype='ag'
end

returninteg_raw(f, inttype, cspec)
end
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
-- integr module follows

localfunctionget_spec(spec, key, mtype)
localdefaults= {eps_abs=1e-6, eps_rel=1e-6, limit=512,
rule='SINGULAR'}
localv=spec[key]
ifvandtype(v) ~=mtypethen
error(string.format('parameter %s should be of type %s', key, mtype))
end
returnvandvordefaults[key]
end

localfunctiontype_check(v, mtype, name)
ifnotvthenerror('parameter ' ..name..' required') end
iftype(v) ~=mtypethen
error('parameter ' ..name..' should be of type ' ..mtype)
end
end

localinteg_raw=integ_module_init ()

functioninteg(spec)
locallimit_min=512

localf=spec.f
type_check (spec.f, 'function', 'f')
ifnotspec.pointsor#spec.points<2then
error('integration points not given')
end
localpts=spec.points
localinttype

localcspec= {}
localinteg_spec_add=
function(key, mtype)
localvalue=get_spec(spec, key, mtype)
cspec[key] =value
end

integ_spec_add ('eps_abs', 'number')
integ_spec_add ('eps_rel', 'number')
integ_spec_add ('limit', 'number')
ifcspec.limit<limit_minthencspec.limit=limit_minend

cspec.a=pts[1]
cspec.b=pts[#pts]

-- infinite intervals
ifpts[1] =='-inf' orpts[#pts] =='+inf' then
ifpts[1] =='-inf' andpts[#pts] =='+inf' then
inttype='agi'
elseifpts[1] =='-inf' then
inttype='agil'
else
inttype='agiu'
end
ifspec.weightthen
localw=spec.weight
ifw.type=='sin' orw.type=='cos' then
inttype=string.gsub(inttype, '^agi', 'awf')
if (inttype=='awf') then
cspec.eps_abs=cspec.eps_abs/2
cspec.a, cspec.b=0, 0
locallint, lerr=integ_raw(f, 'awfl', cspec, w)
localuint, uerr=integ_raw(f, 'awfu', cspec, w)
returnlint+uint, lerr+uerr
end
returninteg_raw(f, inttype, cspec, w)
end
error('cannot calculate indefinite integral with this weight')
end
returninteg_raw(f, inttype, cspec)
end

fori,vinipairs(pts) do
iftype(v) ~='number' thenerror('invalid points') end
end

-- weighted integrals
ifspec.weightthen
localw=spec.weight
if#pts>2then
error('weighted integrals does not accept singular points')
end
ifw.type=='sin' orw.type=='cos' then
type_check (w.omega, 'number', 'omega')
inttype='awo'
elseifw.type=='cauchy' then
type_check (w.singularity, 'number', 'singularity')
inttype='awc'
elseifw.type=='alglog' then
inttype='aws'
else
error('unknown wight type')
end
returninteg_raw(f, inttype, cspec, w)
end

ifnotspec.adaptivethen
if#pts>2then
error('non-adaptive integrals does not accept singular points')
end
inttype='ng'
end

if#pts>2then
integ_spec_add ('points', 'table')
inttype='agp'
end

localrulename=get_spec(spec, 'rule', 'string')
ifrulename=='SINGULAR' then
inttype='ags'
else
cspec.rule=tonumber(string.match(rulename, 'GAUSS(%d+)'))
ifnotcspec.rulethenerror('invalid integration rule') end
inttype='ag'
end

returninteg_raw(f, inttype, cspec)
end