1
;;; Compiled by f2cl version 2.0 beta 2002-05-06
1
;;; Compiled by f2cl version 2.0 beta Date: 2006/01/31 15:11:05
2
;;; Using Lisp CMU Common Lisp Snapshot 2006-01 (19C)
3
4
;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
4
;;; (:coerce-assigns :as-needed) (:array-type ':simple-array)
5
;;; (:coerce-assigns :as-needed) (:array-type ':array)
5
6
;;; (:array-slicing t) (:declare-common nil)
6
7
;;; (:float-format double-float))
11
(let ((x1 (make-array 5 :element-type 'double-float))
12
(w10 (make-array 5 :element-type 'double-float))
13
(x2 (make-array 5 :element-type 'double-float))
14
(w21a (make-array 5 :element-type 'double-float))
15
(w21b (make-array 6 :element-type 'double-float))
16
(x3 (make-array 11 :element-type 'double-float))
17
(w43a (make-array 10 :element-type 'double-float))
18
(w43b (make-array 12 :element-type 'double-float))
19
(x4 (make-array 22 :element-type 'double-float))
20
(w87a (make-array 21 :element-type 'double-float))
21
(w87b (make-array 23 :element-type 'double-float)))
22
(declare (type (simple-array double-float (23)) w87b)
23
(type (simple-array double-float (21)) w87a)
24
(type (simple-array double-float (22)) x4)
25
(type (simple-array double-float (12)) w43b)
26
(type (simple-array double-float (10)) w43a)
27
(type (simple-array double-float (11)) x3)
28
(type (simple-array double-float (6)) w21b)
29
(type (simple-array double-float (5)) w21a x2 w10 x1))
30
(f2cl-lib:fset (f2cl-lib:fref x1 (1) ((1 5))) 0.9739065285171717)
31
(f2cl-lib:fset (f2cl-lib:fref x1 (2) ((1 5))) 0.8650633666889845)
32
(f2cl-lib:fset (f2cl-lib:fref x1 (3) ((1 5))) 0.6794095682990244)
33
(f2cl-lib:fset (f2cl-lib:fref x1 (4) ((1 5))) 0.4333953941292472)
34
(f2cl-lib:fset (f2cl-lib:fref x1 (5) ((1 5))) 0.14887433898163122)
35
(f2cl-lib:fset (f2cl-lib:fref w10 (1) ((1 5))) 0.06667134430868814)
36
(f2cl-lib:fset (f2cl-lib:fref w10 (2) ((1 5))) 0.1494513491505806)
37
(f2cl-lib:fset (f2cl-lib:fref w10 (3) ((1 5))) 0.21908636251598204)
38
(f2cl-lib:fset (f2cl-lib:fref w10 (4) ((1 5))) 0.26926671930999635)
39
(f2cl-lib:fset (f2cl-lib:fref w10 (5) ((1 5))) 0.29552422471475287)
40
(f2cl-lib:fset (f2cl-lib:fref x2 (1) ((1 5))) 0.9956571630258081)
41
(f2cl-lib:fset (f2cl-lib:fref x2 (2) ((1 5))) 0.9301574913557082)
42
(f2cl-lib:fset (f2cl-lib:fref x2 (3) ((1 5))) 0.7808177265864169)
43
(f2cl-lib:fset (f2cl-lib:fref x2 (4) ((1 5))) 0.5627571346686047)
44
(f2cl-lib:fset (f2cl-lib:fref x2 (5) ((1 5))) 0.2943928627014602)
45
(f2cl-lib:fset (f2cl-lib:fref w21a (1) ((1 5))) 0.032558162307964725)
46
(f2cl-lib:fset (f2cl-lib:fref w21a (2) ((1 5))) 0.07503967481091996)
47
(f2cl-lib:fset (f2cl-lib:fref w21a (3) ((1 5))) 0.10938715880229764)
48
(f2cl-lib:fset (f2cl-lib:fref w21a (4) ((1 5))) 0.13470921731147334)
49
(f2cl-lib:fset (f2cl-lib:fref w21a (5) ((1 5))) 0.14773910490133849)
50
(f2cl-lib:fset (f2cl-lib:fref w21b (1) ((1 6))) 0.011694638867371874)
51
(f2cl-lib:fset (f2cl-lib:fref w21b (2) ((1 6))) 0.054755896574351995)
52
(f2cl-lib:fset (f2cl-lib:fref w21b (3) ((1 6))) 0.0931254545836976)
53
(f2cl-lib:fset (f2cl-lib:fref w21b (4) ((1 6))) 0.12349197626206584)
54
(f2cl-lib:fset (f2cl-lib:fref w21b (5) ((1 6))) 0.14277593857706009)
55
(f2cl-lib:fset (f2cl-lib:fref w21b (6) ((1 6))) 0.1494455540029169)
56
(f2cl-lib:fset (f2cl-lib:fref x3 (1) ((1 11))) 0.999333360901932)
57
(f2cl-lib:fset (f2cl-lib:fref x3 (2) ((1 11))) 0.9874334029080889)
58
(f2cl-lib:fset (f2cl-lib:fref x3 (3) ((1 11))) 0.9548079348142663)
59
(f2cl-lib:fset (f2cl-lib:fref x3 (4) ((1 11))) 0.9001486957483283)
60
(f2cl-lib:fset (f2cl-lib:fref x3 (5) ((1 11))) 0.8251983149831141)
61
(f2cl-lib:fset (f2cl-lib:fref x3 (6) ((1 11))) 0.732148388989305)
62
(f2cl-lib:fset (f2cl-lib:fref x3 (7) ((1 11))) 0.6228479705377252)
63
(f2cl-lib:fset (f2cl-lib:fref x3 (8) ((1 11))) 0.4994795740710565)
64
(f2cl-lib:fset (f2cl-lib:fref x3 (9) ((1 11))) 0.36490166134658075)
65
(f2cl-lib:fset (f2cl-lib:fref x3 (10) ((1 11))) 0.2222549197766013)
66
(f2cl-lib:fset (f2cl-lib:fref x3 (11) ((1 11))) 0.07465061746138332)
67
(f2cl-lib:fset (f2cl-lib:fref w43a (1) ((1 10))) 0.016296734289666565)
68
(f2cl-lib:fset (f2cl-lib:fref w43a (2) ((1 10))) 0.0375228761208695)
69
(f2cl-lib:fset (f2cl-lib:fref w43a (3) ((1 10))) 0.05469490205825544)
70
(f2cl-lib:fset (f2cl-lib:fref w43a (4) ((1 10))) 0.06735541460947808)
71
(f2cl-lib:fset (f2cl-lib:fref w43a (5) ((1 10))) 0.07387019963239395)
72
(f2cl-lib:fset (f2cl-lib:fref w43a (6) ((1 10))) 0.005768556059769796)
73
(f2cl-lib:fset (f2cl-lib:fref w43a (7) ((1 10))) 0.027371890593248842)
74
(f2cl-lib:fset (f2cl-lib:fref w43a (8) ((1 10))) 0.04656082691042883)
75
(f2cl-lib:fset (f2cl-lib:fref w43a (9) ((1 10))) 0.06174499520144257)
76
(f2cl-lib:fset (f2cl-lib:fref w43a (10) ((1 10))) 0.07138726726869339)
77
(f2cl-lib:fset (f2cl-lib:fref w43b (1) ((1 12))) 0.001844477640212414)
78
(f2cl-lib:fset (f2cl-lib:fref w43b (2) ((1 12))) 0.010798689585891651)
79
(f2cl-lib:fset (f2cl-lib:fref w43b (3) ((1 12))) 0.021895363867795427)
80
(f2cl-lib:fset (f2cl-lib:fref w43b (4) ((1 12))) 0.032597463975345686)
81
(f2cl-lib:fset (f2cl-lib:fref w43b (5) ((1 12))) 0.04216313793519181)
82
(f2cl-lib:fset (f2cl-lib:fref w43b (6) ((1 12))) 0.050741939600184575)
83
(f2cl-lib:fset (f2cl-lib:fref w43b (7) ((1 12))) 0.05837939554261925)
84
(f2cl-lib:fset (f2cl-lib:fref w43b (8) ((1 12))) 0.06474640495144589)
85
(f2cl-lib:fset (f2cl-lib:fref w43b (9) ((1 12))) 0.06956619791235648)
86
(f2cl-lib:fset (f2cl-lib:fref w43b (10) ((1 12))) 0.07282444147183322)
87
(f2cl-lib:fset (f2cl-lib:fref w43b (11) ((1 12))) 0.07450775101417512)
88
(f2cl-lib:fset (f2cl-lib:fref w43b (12) ((1 12))) 0.07472214751740301)
89
(f2cl-lib:fset (f2cl-lib:fref x4 (1) ((1 22))) 0.9999029772627293)
90
(f2cl-lib:fset (f2cl-lib:fref x4 (2) ((1 22))) 0.9979898959866788)
91
(f2cl-lib:fset (f2cl-lib:fref x4 (3) ((1 22))) 0.9921754978606873)
92
(f2cl-lib:fset (f2cl-lib:fref x4 (4) ((1 22))) 0.9813581635727128)
93
(f2cl-lib:fset (f2cl-lib:fref x4 (5) ((1 22))) 0.9650576238583847)
94
(f2cl-lib:fset (f2cl-lib:fref x4 (6) ((1 22))) 0.9431676131336706)
95
(f2cl-lib:fset (f2cl-lib:fref x4 (7) ((1 22))) 0.9158064146855072)
96
(f2cl-lib:fset (f2cl-lib:fref x4 (8) ((1 22))) 0.8832216577713164)
97
(f2cl-lib:fset (f2cl-lib:fref x4 (9) ((1 22))) 0.8457107484624157)
98
(f2cl-lib:fset (f2cl-lib:fref x4 (10) ((1 22))) 0.8035576580352309)
99
(f2cl-lib:fset (f2cl-lib:fref x4 (11) ((1 22))) 0.7570057306854956)
100
(f2cl-lib:fset (f2cl-lib:fref x4 (12) ((1 22))) 0.7062732097873218)
101
(f2cl-lib:fset (f2cl-lib:fref x4 (13) ((1 22))) 0.6515894665011779)
102
(f2cl-lib:fset (f2cl-lib:fref x4 (14) ((1 22))) 0.5932233740579611)
103
(f2cl-lib:fset (f2cl-lib:fref x4 (15) ((1 22))) 0.531493605970832)
104
(f2cl-lib:fset (f2cl-lib:fref x4 (16) ((1 22))) 0.46676362304202285)
105
(f2cl-lib:fset (f2cl-lib:fref x4 (17) ((1 22))) 0.3994248478592188)
106
(f2cl-lib:fset (f2cl-lib:fref x4 (18) ((1 22))) 0.3298748771061883)
107
(f2cl-lib:fset (f2cl-lib:fref x4 (19) ((1 22))) 0.25850355920216156)
108
(f2cl-lib:fset (f2cl-lib:fref x4 (20) ((1 22))) 0.18569539656834666)
109
(f2cl-lib:fset (f2cl-lib:fref x4 (21) ((1 22))) 0.11184221317990747)
110
(f2cl-lib:fset (f2cl-lib:fref x4 (22) ((1 22))) 0.03735212339461987)
111
(f2cl-lib:fset (f2cl-lib:fref w87a (1) ((1 21))) 0.008148377384149173)
112
(f2cl-lib:fset (f2cl-lib:fref w87a (2) ((1 21))) 0.018761438201562824)
113
(f2cl-lib:fset (f2cl-lib:fref w87a (3) ((1 21))) 0.027347451050052287)
114
(f2cl-lib:fset (f2cl-lib:fref w87a (4) ((1 21))) 0.03367770731163793)
115
(f2cl-lib:fset (f2cl-lib:fref w87a (5) ((1 21))) 0.036935099820427905)
116
(f2cl-lib:fset (f2cl-lib:fref w87a (6) ((1 21))) 0.0028848724302115306)
117
(f2cl-lib:fset (f2cl-lib:fref w87a (7) ((1 21))) 0.013685946022712702)
118
(f2cl-lib:fset (f2cl-lib:fref w87a (8) ((1 21))) 0.02328041350288831)
119
(f2cl-lib:fset (f2cl-lib:fref w87a (9) ((1 21))) 0.03087249761171336)
120
(f2cl-lib:fset (f2cl-lib:fref w87a (10) ((1 21))) 0.03569363363941877)
121
(f2cl-lib:fset (f2cl-lib:fref w87a (11) ((1 21))) 9.152833452022413e-4)
122
(f2cl-lib:fset (f2cl-lib:fref w87a (12) ((1 21))) 0.005399280219300471)
123
(f2cl-lib:fset (f2cl-lib:fref w87a (13) ((1 21))) 0.01094767960111893)
124
(f2cl-lib:fset (f2cl-lib:fref w87a (14) ((1 21))) 0.016298731696787336)
125
(f2cl-lib:fset (f2cl-lib:fref w87a (15) ((1 21))) 0.021081568889203834)
126
(f2cl-lib:fset (f2cl-lib:fref w87a (16) ((1 21))) 0.025370969769253827)
127
(f2cl-lib:fset (f2cl-lib:fref w87a (17) ((1 21))) 0.029189697756475754)
128
(f2cl-lib:fset (f2cl-lib:fref w87a (18) ((1 21))) 0.03237320246720279)
129
(f2cl-lib:fset (f2cl-lib:fref w87a (19) ((1 21))) 0.034783098950365146)
130
(f2cl-lib:fset (f2cl-lib:fref w87a (20) ((1 21))) 0.03641222073135179)
131
(f2cl-lib:fset (f2cl-lib:fref w87a (21) ((1 21))) 0.037253875503047706)
132
(f2cl-lib:fset (f2cl-lib:fref w87b (1) ((1 23))) 2.7414556376207233e-4)
133
(f2cl-lib:fset (f2cl-lib:fref w87b (2) ((1 23))) 0.0018071241550579428)
134
(f2cl-lib:fset (f2cl-lib:fref w87b (3) ((1 23))) 0.0040968692827591646)
135
(f2cl-lib:fset (f2cl-lib:fref w87b (4) ((1 23))) 0.006758290051847379)
136
(f2cl-lib:fset (f2cl-lib:fref w87b (5) ((1 23))) 0.009549957672201646)
137
(f2cl-lib:fset (f2cl-lib:fref w87b (6) ((1 23))) 0.012329447652244854)
138
(f2cl-lib:fset (f2cl-lib:fref w87b (7) ((1 23))) 0.015010447346388952)
139
(f2cl-lib:fset (f2cl-lib:fref w87b (8) ((1 23))) 0.01754896798624319)
140
(f2cl-lib:fset (f2cl-lib:fref w87b (9) ((1 23))) 0.019938037786440887)
141
(f2cl-lib:fset (f2cl-lib:fref w87b (10) ((1 23))) 0.022194935961012286)
142
(f2cl-lib:fset (f2cl-lib:fref w87b (11) ((1 23))) 0.024339147126000805)
143
(f2cl-lib:fset (f2cl-lib:fref w87b (12) ((1 23))) 0.026374505414839208)
144
(f2cl-lib:fset (f2cl-lib:fref w87b (13) ((1 23))) 0.0282869107887712)
145
(f2cl-lib:fset (f2cl-lib:fref w87b (14) ((1 23))) 0.030052581128092695)
146
(f2cl-lib:fset (f2cl-lib:fref w87b (15) ((1 23))) 0.03164675137143993)
147
(f2cl-lib:fset (f2cl-lib:fref w87b (16) ((1 23))) 0.033050413419978504)
148
(f2cl-lib:fset (f2cl-lib:fref w87b (17) ((1 23))) 0.034255099704226064)
149
(f2cl-lib:fset (f2cl-lib:fref w87b (18) ((1 23))) 0.03526241266015668)
150
(f2cl-lib:fset (f2cl-lib:fref w87b (19) ((1 23))) 0.0360769896228887)
151
(f2cl-lib:fset (f2cl-lib:fref w87b (20) ((1 23))) 0.03669860449845609)
152
(f2cl-lib:fset (f2cl-lib:fref w87b (21) ((1 23))) 0.037120549269832576)
153
(f2cl-lib:fset (f2cl-lib:fref w87b (22) ((1 23))) 0.03733422875193504)
154
(f2cl-lib:fset (f2cl-lib:fref w87b (23) ((1 23))) 0.037361073762679026)
14
:element-type 'double-float
15
:initial-contents '(0.9739065285171717 0.8650633666889845
16
0.6794095682990244 0.4333953941292472
17
0.14887433898163122)))
20
:element-type 'double-float
21
:initial-contents '(0.06667134430868814 0.1494513491505806
22
0.21908636251598204 0.26926671930999635
23
0.29552422471475287)))
26
:element-type 'double-float
27
:initial-contents '(0.9956571630258081 0.9301574913557082
28
0.7808177265864169 0.5627571346686047
32
:element-type 'double-float
33
:initial-contents '(0.032558162307964725 0.07503967481091996
34
0.10938715880229764 0.13470921731147334
35
0.14773910490133849)))
38
:element-type 'double-float
39
:initial-contents '(0.011694638867371874
40
0.054755896574351995 0.0931254545836976
41
0.12349197626206584 0.14277593857706009
45
:element-type 'double-float
46
:initial-contents '(0.999333360901932 0.9874334029080889
47
0.9548079348142663 0.9001486957483283
48
0.8251983149831141 0.732148388989305
49
0.6228479705377252 0.4994795740710565
50
0.36490166134658075 0.2222549197766013
51
0.07465061746138332)))
54
:element-type 'double-float
55
:initial-contents '(0.016296734289666565 0.0375228761208695
56
0.05469490205825544 0.06735541460947808
57
0.07387019963239395 0.005768556059769796
58
0.027371890593248842 0.04656082691042883
60
0.07138726726869339)))
63
:element-type 'double-float
64
:initial-contents '(0.001844477640212414
67
0.032597463975345686 0.04216313793519181
68
0.050741939600184575 0.05837939554261925
69
0.06474640495144589 0.06956619791235648
70
0.07282444147183322 0.07450775101417512
71
0.07472214751740301)))
74
:element-type 'double-float
75
:initial-contents '(0.9999029772627293 0.9979898959866788
76
0.9921754978606873 0.9813581635727128
77
0.9650576238583847 0.9431676131336706
78
0.9158064146855072 0.8832216577713164
79
0.8457107484624157 0.8035576580352309
80
0.7570057306854956 0.7062732097873218
81
0.6515894665011779 0.5932233740579611
82
0.531493605970832 0.46676362304202285
83
0.3994248478592188 0.3298748771061883
84
0.25850355920216156 0.18569539656834666
86
0.03735212339461987)))
89
:element-type 'double-float
90
:initial-contents '(0.008148377384149173
92
0.027347451050052287 0.03367770731163793
95
0.013685946022712702 0.02328041350288831
96
0.03087249761171336 0.03569363363941877
98
0.005399280219300471 0.01094767960111893
102
0.029189697756475754 0.03237320246720279
103
0.034783098950365146 0.03641222073135179
104
0.037253875503047706)))
107
:element-type 'double-float
108
:initial-contents '(2.7414556376207234e-4
109
0.0018071241550579428
110
0.0040968692827591646
114
0.015010447346388952 0.01754896798624319
118
0.026374505414839208 0.0282869107887712
119
0.030052581128092695 0.03164675137143993
121
0.034255099704226064 0.03526241266015668
122
0.0360769896228887 0.03669860449845609
123
0.037120549269832576 0.03733422875193504
124
0.037361073762679026))))
125
(declare (type (array double-float (23)) w87b)
126
(type (array double-float (21)) w87a)
127
(type (array double-float (22)) x4)
128
(type (array double-float (12)) w43b)
129
(type (array double-float (10)) w43a)
130
(type (array double-float (11)) x3)
131
(type (array double-float (6)) w21b)
132
(type (array double-float (5)) w21a x2 w10 x1))
155
133
(defun dqng (f a b epsabs epsrel result abserr neval ier)
156
134
(declare (type f2cl-lib:integer4 ier neval)
157
(type double-float abserr result epsrel epsabs b a)
158
(type (function (double-float) (values double-float &rest t)) f))
159
(prog ((fv1 (make-array 5 :element-type 'double-float))
160
(fv2 (make-array 5 :element-type 'double-float))
161
(fv3 (make-array 5 :element-type 'double-float))
162
(fv4 (make-array 5 :element-type 'double-float))
163
(savfun (make-array 21 :element-type 'double-float)) (ipx 0) (k 0)
164
(l 0) (absc 0.0) (centr 0.0) (dhlgth 0.0) (epmach 0.0) (fcentr 0.0)
165
(fval 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (res10 0.0)
166
(res21 0.0) (res43 0.0) (res87 0.0) (resabs 0.0) (resasc 0.0)
167
(reskh 0.0) (uflow 0.0) (abs$ 0.0f0))
168
(declare (type single-float abs$)
169
(type (simple-array double-float (21)) savfun)
170
(type (simple-array double-float (5)) fv4 fv3 fv2 fv1)
171
(type double-float uflow reskh resasc resabs res87 res43 res21 res10
172
hlgth fval2 fval1 fval fcentr epmach dhlgth centr absc)
173
(type f2cl-lib:integer4 l k ipx))
174
(setf epmach (f2cl-lib:d1mach 4))
175
(setf uflow (f2cl-lib:d1mach 1))
180
(if (and (<= epsabs 0.0) (< epsrel (max (* 50.0 epmach) 5.0e-29)))
182
(setf hlgth (* 0.5 (- b a)))
183
(setf dhlgth (coerce (abs hlgth) 'double-float))
184
(setf centr (* 0.5 (+ b a)))
190
(when var-0 (setf centr var-0))
194
(f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
197
(f2cl-lib:computed-goto (label5 label25 label45) l)
200
(setf res21 (* (f2cl-lib:fref w21b (6) ((1 6))) fcentr))
201
(setf resabs (* (f2cl-lib:fref w21b (6) ((1 6))) (abs fcentr)))
202
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
205
(setf absc (* hlgth (f2cl-lib:fref x1 (k) ((1 5)))))
206
(setf fval1 (funcall f (+ centr absc)))
207
(setf fval2 (funcall f (- centr absc)))
208
(setf fval (+ fval1 fval2))
209
(setf res10 (+ res10 (* (f2cl-lib:fref w10 (k) ((1 5))) fval)))
210
(setf res21 (+ res21 (* (f2cl-lib:fref w21a (k) ((1 5))) fval)))
213
(* (f2cl-lib:fref w21a (k) ((1 5)))
214
(+ (abs fval1) (abs fval2)))))
215
(f2cl-lib:fset (f2cl-lib:fref savfun (k) ((1 21))) fval)
216
(f2cl-lib:fset (f2cl-lib:fref fv1 (k) ((1 5))) fval1)
217
(f2cl-lib:fset (f2cl-lib:fref fv2 (k) ((1 5))) fval2)
220
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
223
(setf ipx (f2cl-lib:int-add ipx 1))
224
(setf absc (* hlgth (f2cl-lib:fref x2 (k) ((1 5)))))
225
(setf fval1 (funcall f (+ centr absc)))
226
(setf fval2 (funcall f (- centr absc)))
227
(setf fval (+ fval1 fval2))
228
(setf res21 (+ res21 (* (f2cl-lib:fref w21b (k) ((1 6))) fval)))
231
(* (f2cl-lib:fref w21b (k) ((1 6)))
232
(+ (abs fval1) (abs fval2)))))
233
(f2cl-lib:fset (f2cl-lib:fref savfun (ipx) ((1 21))) fval)
234
(f2cl-lib:fset (f2cl-lib:fref fv3 (k) ((1 5))) fval1)
235
(f2cl-lib:fset (f2cl-lib:fref fv4 (k) ((1 5))) fval2)
237
(setf result (* res21 hlgth))
238
(setf resabs (* resabs dhlgth))
239
(setf reskh (* 0.5 res21))
241
(* (f2cl-lib:fref w21b (6) ((1 6))) (abs (- fcentr reskh))))
242
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
247
(* (f2cl-lib:fref w21a (k) ((1 5)))
248
(+ (abs (- (f2cl-lib:fref fv1 (k) ((1 5))) reskh))
135
(type double-float abserr result epsrel epsabs b a))
136
(f2cl-lib:with-multi-array-data
138
(prog ((fv1 (make-array 5 :element-type 'double-float))
139
(fv2 (make-array 5 :element-type 'double-float))
140
(fv3 (make-array 5 :element-type 'double-float))
141
(fv4 (make-array 5 :element-type 'double-float))
142
(savfun (make-array 21 :element-type 'double-float)) (ipx 0) (k 0)
143
(l 0) (absc 0.0) (centr 0.0) (dhlgth 0.0) (epmach 0.0)
144
(fcentr 0.0) (fval 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0)
145
(res10 0.0) (res21 0.0) (res43 0.0) (res87 0.0) (resabs 0.0)
146
(resasc 0.0) (reskh 0.0) (uflow 0.0) (abs$ 0.0f0))
147
(declare (type single-float abs$)
148
(type (array double-float (21)) savfun)
149
(type (array double-float (5)) fv4 fv3 fv2 fv1)
150
(type double-float uflow reskh resasc resabs res87 res43 res21
151
res10 hlgth fval2 fval1 fval fcentr epmach
153
(type f2cl-lib:integer4 l k ipx))
154
(setf epmach (f2cl-lib:d1mach 4))
155
(setf uflow (f2cl-lib:d1mach 1))
160
(if (and (<= epsabs 0.0) (< epsrel (max (* 50.0 epmach) 5.e-29)))
162
(setf hlgth (* 0.5 (- b a)))
163
(setf dhlgth (coerce (abs hlgth) 'double-float))
164
(setf centr (* 0.5 (+ b a)))
166
(multiple-value-bind (ret-val var-0)
174
(f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
177
(f2cl-lib:computed-goto (label5 label25 label45) l)
180
(setf res21 (* (f2cl-lib:fref w21b (6) ((1 6))) fcentr))
181
(setf resabs (* (f2cl-lib:fref w21b (6) ((1 6))) (abs fcentr)))
182
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
185
(setf absc (* hlgth (f2cl-lib:fref x1 (k) ((1 5)))))
186
(setf fval1 (funcall f (+ centr absc)))
187
(setf fval2 (funcall f (- centr absc)))
188
(setf fval (+ fval1 fval2))
189
(setf res10 (+ res10 (* (f2cl-lib:fref w10 (k) ((1 5))) fval)))
191
(+ res21 (* (f2cl-lib:fref w21a (k) ((1 5))) fval)))
194
(* (f2cl-lib:fref w21a (k) ((1 5)))
195
(+ (abs fval1) (abs fval2)))))
196
(f2cl-lib:fset (f2cl-lib:fref savfun (k) ((1 21))) fval)
197
(f2cl-lib:fset (f2cl-lib:fref fv1 (k) ((1 5))) fval1)
198
(f2cl-lib:fset (f2cl-lib:fref fv2 (k) ((1 5))) fval2)
201
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
204
(setf ipx (f2cl-lib:int-add ipx 1))
205
(setf absc (* hlgth (f2cl-lib:fref x2 (k) ((1 5)))))
206
(setf fval1 (funcall f (+ centr absc)))
207
(setf fval2 (funcall f (- centr absc)))
208
(setf fval (+ fval1 fval2))
210
(+ res21 (* (f2cl-lib:fref w21b (k) ((1 6))) fval)))
213
(* (f2cl-lib:fref w21b (k) ((1 6)))
214
(+ (abs fval1) (abs fval2)))))
215
(f2cl-lib:fset (f2cl-lib:fref savfun (ipx) ((1 21))) fval)
216
(f2cl-lib:fset (f2cl-lib:fref fv3 (k) ((1 5))) fval1)
217
(f2cl-lib:fset (f2cl-lib:fref fv4 (k) ((1 5))) fval2)
219
(setf result (* res21 hlgth))
220
(setf resabs (* resabs dhlgth))
221
(setf reskh (* 0.5 res21))
223
(* (f2cl-lib:fref w21b (6) ((1 6)))
224
(abs (- fcentr reskh))))
225
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
230
(* (f2cl-lib:fref w21a (k) ((1 5)))
232
(abs (- (f2cl-lib:fref fv1 (k) ((1 5))) reskh))
250
234
(- (f2cl-lib:fref fv2 (k) ((1 5))) reskh))))
251
(* (f2cl-lib:fref w21b (k) ((1 6)))
252
(+ (abs (- (f2cl-lib:fref fv3 (k) ((1 5))) reskh))
235
(* (f2cl-lib:fref w21b (k) ((1 6)))
237
(abs (- (f2cl-lib:fref fv3 (k) ((1 5))) reskh))
254
239
(- (f2cl-lib:fref fv4 (k) ((1 5))) reskh))))))
256
(setf abserr (coerce (abs (* (- res21 res10) hlgth)) 'double-float))
257
(setf resasc (* resasc dhlgth))
260
(setf res43 (* (f2cl-lib:fref w43b (12) ((1 12))) fcentr))
262
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
267
(* (f2cl-lib:fref savfun (k) ((1 21)))
268
(f2cl-lib:fref w43a (k) ((1 10))))))
270
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
273
(setf ipx (f2cl-lib:int-add ipx 1))
274
(setf absc (* hlgth (f2cl-lib:fref x3 (k) ((1 11)))))
276
(+ (funcall f (+ absc centr))
277
(funcall f (- centr absc))))
278
(setf res43 (+ res43 (* fval (f2cl-lib:fref w43b (k) ((1 12))))))
279
(f2cl-lib:fset (f2cl-lib:fref savfun (ipx) ((1 21))) fval)
281
(setf result (* res43 hlgth))
282
(setf abserr (coerce (abs (* (- res43 res21) hlgth)) 'double-float))
285
(setf res87 (* (f2cl-lib:fref w87b (23) ((1 23))) fcentr))
287
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
292
(* (f2cl-lib:fref savfun (k) ((1 21)))
293
(f2cl-lib:fref w87a (k) ((1 21))))))
295
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
298
(setf absc (* hlgth (f2cl-lib:fref x4 (k) ((1 22)))))
301
(* (f2cl-lib:fref w87b (k) ((1 23)))
302
(+ (funcall f (+ absc centr))
303
(funcall f (- centr absc))))))
305
(setf result (* res87 hlgth))
306
(setf abserr (coerce (abs (* (- res87 res43) hlgth)) 'double-float))
308
(if (and (/= resasc 0.0) (/= abserr 0.0))
311
(min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
312
(if (> resabs (/ uflow (* 50.0 epmach)))
313
(setf abserr (max (* epmach 50.0 resabs) abserr)))
314
(if (<= abserr (max epsabs (* epsrel (abs result)))) (setf ier 0))
315
(if (= ier 0) (go label999))
319
(var-0 var-1 var-2 var-3 var-4)
320
(xermsg "SLATEC" "DQNG" "ABNORMAL RETURN" ier 0)
321
(declare (ignore var-0 var-1 var-2 var-4))
322
(when var-3 (setf ier var-3)))
326
(return (values nil nil nil nil nil result abserr neval ier)))))
242
(coerce (abs (* (- res21 res10) hlgth)) 'double-float))
243
(setf resasc (* resasc dhlgth))
246
(setf res43 (* (f2cl-lib:fref w43b (12) ((1 12))) fcentr))
248
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
253
(* (f2cl-lib:fref savfun (k) ((1 21)))
254
(f2cl-lib:fref w43a (k) ((1 10))))))
256
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
259
(setf ipx (f2cl-lib:int-add ipx 1))
260
(setf absc (* hlgth (f2cl-lib:fref x3 (k) ((1 11)))))
262
(+ (funcall f (+ absc centr))
263
(funcall f (- centr absc))))
265
(+ res43 (* fval (f2cl-lib:fref w43b (k) ((1 12))))))
266
(f2cl-lib:fset (f2cl-lib:fref savfun (ipx) ((1 21))) fval)
268
(setf result (* res43 hlgth))
270
(coerce (abs (* (- res43 res21) hlgth)) 'double-float))
273
(setf res87 (* (f2cl-lib:fref w87b (23) ((1 23))) fcentr))
275
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
280
(* (f2cl-lib:fref savfun (k) ((1 21)))
281
(f2cl-lib:fref w87a (k) ((1 21))))))
283
(f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
286
(setf absc (* hlgth (f2cl-lib:fref x4 (k) ((1 22)))))
289
(* (f2cl-lib:fref w87b (k) ((1 23)))
290
(+ (funcall f (+ absc centr))
291
(funcall f (- centr absc))))))
293
(setf result (* res87 hlgth))
295
(coerce (abs (* (- res87 res43) hlgth)) 'double-float))
297
(if (and (/= resasc 0.0) (/= abserr 0.0))
300
(min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
301
(if (> resabs (/ uflow (* 50.0 epmach)))
302
(setf abserr (max (* epmach 50.0 resabs) abserr)))
303
(if (<= abserr (max epsabs (* epsrel (abs result)))) (setf ier 0))
304
(if (= ier 0) (go label999))
307
(xermsg "SLATEC" "DQNG" "ABNORMAL RETURN" ier 0)
311
(return (values nil nil nil nil nil result abserr neval ier))))))