~ubuntu-branches/debian/squeeze/maxima/squeeze

« back to all changes in this revision

Viewing changes to src/numerical/slatec/dqng.lisp

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2006-10-18 14:52:42 UTC
  • mto: (1.1.5 upstream)
  • mto: This revision was merged to the branch mainline in revision 4.
  • Revision ID: james.westby@ubuntu.com-20061018145242-vzyrm5hmxr8kiosf
ImportĀ upstreamĀ versionĀ 5.10.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
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)
2
3
;;; 
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))
7
8
 
8
 
(in-package "SLATEC")
9
 
 
10
 
 
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)
 
9
(in-package :slatec)
 
10
 
 
11
 
 
12
(let ((x1
 
13
       (make-array 5
 
14
                   :element-type 'double-float
 
15
                   :initial-contents '(0.9739065285171717 0.8650633666889845
 
16
                                       0.6794095682990244 0.4333953941292472
 
17
                                       0.14887433898163122)))
 
18
      (w10
 
19
       (make-array 5
 
20
                   :element-type 'double-float
 
21
                   :initial-contents '(0.06667134430868814 0.1494513491505806
 
22
                                       0.21908636251598204 0.26926671930999635
 
23
                                       0.29552422471475287)))
 
24
      (x2
 
25
       (make-array 5
 
26
                   :element-type 'double-float
 
27
                   :initial-contents '(0.9956571630258081 0.9301574913557082
 
28
                                       0.7808177265864169 0.5627571346686047
 
29
                                       0.2943928627014602)))
 
30
      (w21a
 
31
       (make-array 5
 
32
                   :element-type 'double-float
 
33
                   :initial-contents '(0.032558162307964725 0.07503967481091996
 
34
                                       0.10938715880229764 0.13470921731147334
 
35
                                       0.14773910490133849)))
 
36
      (w21b
 
37
       (make-array 6
 
38
                   :element-type 'double-float
 
39
                   :initial-contents '(0.011694638867371874
 
40
                                       0.054755896574351995 0.0931254545836976
 
41
                                       0.12349197626206584 0.14277593857706009
 
42
                                       0.1494455540029169)))
 
43
      (x3
 
44
       (make-array 11
 
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)))
 
52
      (w43a
 
53
       (make-array 10
 
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
 
59
                                       0.06174499520144257
 
60
                                       0.07138726726869339)))
 
61
      (w43b
 
62
       (make-array 12
 
63
                   :element-type 'double-float
 
64
                   :initial-contents '(0.001844477640212414
 
65
                                       0.010798689585891651
 
66
                                       0.021895363867795427
 
67
                                       0.032597463975345686 0.04216313793519181
 
68
                                       0.050741939600184575 0.05837939554261925
 
69
                                       0.06474640495144589 0.06956619791235648
 
70
                                       0.07282444147183322 0.07450775101417512
 
71
                                       0.07472214751740301)))
 
72
      (x4
 
73
       (make-array 22
 
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
 
85
                                       0.11184221317990747
 
86
                                       0.03735212339461987)))
 
87
      (w87a
 
88
       (make-array 21
 
89
                   :element-type 'double-float
 
90
                   :initial-contents '(0.008148377384149173
 
91
                                       0.018761438201562824
 
92
                                       0.027347451050052287 0.03367770731163793
 
93
                                       0.036935099820427905
 
94
                                       0.0028848724302115306
 
95
                                       0.013685946022712702 0.02328041350288831
 
96
                                       0.03087249761171336 0.03569363363941877
 
97
                                       9.152833452022414e-4
 
98
                                       0.005399280219300471 0.01094767960111893
 
99
                                       0.016298731696787336
 
100
                                       0.021081568889203834
 
101
                                       0.025370969769253827
 
102
                                       0.029189697756475754 0.03237320246720279
 
103
                                       0.034783098950365146 0.03641222073135179
 
104
                                       0.037253875503047706)))
 
105
      (w87b
 
106
       (make-array 23
 
107
                   :element-type 'double-float
 
108
                   :initial-contents '(2.7414556376207234e-4
 
109
                                       0.0018071241550579428
 
110
                                       0.0040968692827591646
 
111
                                       0.006758290051847379
 
112
                                       0.009549957672201646
 
113
                                       0.012329447652244854
 
114
                                       0.015010447346388952 0.01754896798624319
 
115
                                       0.019938037786440887
 
116
                                       0.022194935961012286
 
117
                                       0.024339147126000805
 
118
                                       0.026374505414839208 0.0282869107887712
 
119
                                       0.030052581128092695 0.03164675137143993
 
120
                                       0.033050413419978504
 
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))
176
 
      (setf result 0.0)
177
 
      (setf abserr 0.0)
178
 
      (setf neval 0)
179
 
      (setf ier 6)
180
 
      (if (and (<= epsabs 0.0) (< epsrel (max (* 50.0 epmach) 5.0e-29)))
181
 
          (go label80))
182
 
      (setf hlgth (* 0.5 (- b a)))
183
 
      (setf dhlgth (coerce (abs hlgth) 'double-float))
184
 
      (setf centr (* 0.5 (+ b a)))
185
 
      (setf fcentr
186
 
              (multiple-value-bind
187
 
                  (ret-val var-0)
188
 
                  (funcall f centr)
189
 
                (declare (ignore))
190
 
                (when var-0 (setf centr var-0))
191
 
                ret-val))
192
 
      (setf neval 21)
193
 
      (setf ier 1)
194
 
      (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
195
 
                    ((> l 3) nil)
196
 
        (tagbody
197
 
          (f2cl-lib:computed-goto (label5 label25 label45) l)
198
 
         label5
199
 
          (setf res10 0.0)
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))
203
 
                        ((> k 5) nil)
204
 
            (tagbody
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)))
211
 
              (setf resabs
212
 
                      (+ resabs
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)
218
 
             label10))
219
 
          (setf ipx 5)
220
 
          (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
221
 
                        ((> k 5) nil)
222
 
            (tagbody
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)))
229
 
              (setf resabs
230
 
                      (+ resabs
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)
236
 
             label15))
237
 
          (setf result (* res21 hlgth))
238
 
          (setf resabs (* resabs dhlgth))
239
 
          (setf reskh (* 0.5 res21))
240
 
          (setf resasc
241
 
                  (* (f2cl-lib:fref w21b (6) ((1 6))) (abs (- fcentr reskh))))
242
 
          (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
243
 
                        ((> k 5) nil)
244
 
            (tagbody
245
 
              (setf resasc
246
 
                      (+ resasc
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
 
137
        nil
 
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
 
152
                                    dhlgth centr absc)
 
153
                 (type f2cl-lib:integer4 l k ipx))
 
154
        (setf epmach (f2cl-lib:d1mach 4))
 
155
        (setf uflow (f2cl-lib:d1mach 1))
 
156
        (setf result 0.0)
 
157
        (setf abserr 0.0)
 
158
        (setf neval 0)
 
159
        (setf ier 6)
 
160
        (if (and (<= epsabs 0.0) (< epsrel (max (* 50.0 epmach) 5.e-29)))
 
161
            (go label80))
 
162
        (setf hlgth (* 0.5 (- b a)))
 
163
        (setf dhlgth (coerce (abs hlgth) 'double-float))
 
164
        (setf centr (* 0.5 (+ b a)))
 
165
        (setf fcentr
 
166
                (multiple-value-bind (ret-val var-0)
 
167
                    (funcall f centr)
 
168
                  (declare (ignore))
 
169
                  (when var-0
 
170
                    (setf centr var-0))
 
171
                  ret-val))
 
172
        (setf neval 21)
 
173
        (setf ier 1)
 
174
        (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 
175
                      ((> l 3) nil)
 
176
          (tagbody
 
177
            (f2cl-lib:computed-goto (label5 label25 label45) l)
 
178
           label5
 
179
            (setf res10 0.0)
 
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))
 
183
                          ((> k 5) nil)
 
184
              (tagbody
 
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)))
 
190
                (setf res21
 
191
                        (+ res21 (* (f2cl-lib:fref w21a (k) ((1 5))) fval)))
 
192
                (setf resabs
 
193
                        (+ resabs
 
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)
 
199
               label10))
 
200
            (setf ipx 5)
 
201
            (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 
202
                          ((> k 5) nil)
 
203
              (tagbody
 
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))
 
209
                (setf res21
 
210
                        (+ res21 (* (f2cl-lib:fref w21b (k) ((1 6))) fval)))
 
211
                (setf resabs
 
212
                        (+ resabs
 
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)
 
218
               label15))
 
219
            (setf result (* res21 hlgth))
 
220
            (setf resabs (* resabs dhlgth))
 
221
            (setf reskh (* 0.5 res21))
 
222
            (setf resasc
 
223
                    (* (f2cl-lib:fref w21b (6) ((1 6)))
 
224
                       (abs (- fcentr reskh))))
 
225
            (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 
226
                          ((> k 5) nil)
 
227
              (tagbody
 
228
                (setf resasc
 
229
                        (+ resasc
 
230
                           (* (f2cl-lib:fref w21a (k) ((1 5)))
 
231
                              (+
 
232
                               (abs (- (f2cl-lib:fref fv1 (k) ((1 5))) reskh))
249
233
                               (abs
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)))
 
236
                              (+
 
237
                               (abs (- (f2cl-lib:fref fv3 (k) ((1 5))) reskh))
253
238
                               (abs
254
239
                                (- (f2cl-lib:fref fv4 (k) ((1 5))) reskh))))))
255
 
             label20))
256
 
          (setf abserr (coerce (abs (* (- res21 res10) hlgth)) 'double-float))
257
 
          (setf resasc (* resasc dhlgth))
258
 
          (go label65)
259
 
         label25
260
 
          (setf res43 (* (f2cl-lib:fref w43b (12) ((1 12))) fcentr))
261
 
          (setf neval 43)
262
 
          (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
263
 
                        ((> k 10) nil)
264
 
            (tagbody
265
 
              (setf res43
266
 
                      (+ res43
267
 
                         (* (f2cl-lib:fref savfun (k) ((1 21)))
268
 
                            (f2cl-lib:fref w43a (k) ((1 10))))))
269
 
             label30))
270
 
          (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
271
 
                        ((> k 11) nil)
272
 
            (tagbody
273
 
              (setf ipx (f2cl-lib:int-add ipx 1))
274
 
              (setf absc (* hlgth (f2cl-lib:fref x3 (k) ((1 11)))))
275
 
              (setf fval
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)
280
 
             label40))
281
 
          (setf result (* res43 hlgth))
282
 
          (setf abserr (coerce (abs (* (- res43 res21) hlgth)) 'double-float))
283
 
          (go label65)
284
 
         label45
285
 
          (setf res87 (* (f2cl-lib:fref w87b (23) ((1 23))) fcentr))
286
 
          (setf neval 87)
287
 
          (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
288
 
                        ((> k 21) nil)
289
 
            (tagbody
290
 
              (setf res87
291
 
                      (+ res87
292
 
                         (* (f2cl-lib:fref savfun (k) ((1 21)))
293
 
                            (f2cl-lib:fref w87a (k) ((1 21))))))
294
 
             label50))
295
 
          (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
296
 
                        ((> k 22) nil)
297
 
            (tagbody
298
 
              (setf absc (* hlgth (f2cl-lib:fref x4 (k) ((1 22)))))
299
 
              (setf res87
300
 
                      (+ res87
301
 
                         (* (f2cl-lib:fref w87b (k) ((1 23)))
302
 
                            (+ (funcall f (+ absc centr))
303
 
                               (funcall f (- centr absc))))))
304
 
             label60))
305
 
          (setf result (* res87 hlgth))
306
 
          (setf abserr (coerce (abs (* (- res87 res43) hlgth)) 'double-float))
307
 
         label65
308
 
          (if (and (/= resasc 0.0) (/= abserr 0.0))
309
 
              (setf abserr
310
 
                      (* resasc
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))
316
 
         label70))
317
 
     label80
318
 
      (multiple-value-bind
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)))
323
 
     label999
324
 
      (go end_label)
325
 
     end_label
326
 
      (return (values nil nil nil nil nil result abserr neval ier)))))
 
240
               label20))
 
241
            (setf abserr
 
242
                    (coerce (abs (* (- res21 res10) hlgth)) 'double-float))
 
243
            (setf resasc (* resasc dhlgth))
 
244
            (go label65)
 
245
           label25
 
246
            (setf res43 (* (f2cl-lib:fref w43b (12) ((1 12))) fcentr))
 
247
            (setf neval 43)
 
248
            (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 
249
                          ((> k 10) nil)
 
250
              (tagbody
 
251
                (setf res43
 
252
                        (+ res43
 
253
                           (* (f2cl-lib:fref savfun (k) ((1 21)))
 
254
                              (f2cl-lib:fref w43a (k) ((1 10))))))
 
255
               label30))
 
256
            (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 
257
                          ((> k 11) nil)
 
258
              (tagbody
 
259
                (setf ipx (f2cl-lib:int-add ipx 1))
 
260
                (setf absc (* hlgth (f2cl-lib:fref x3 (k) ((1 11)))))
 
261
                (setf fval
 
262
                        (+ (funcall f (+ absc centr))
 
263
                           (funcall f (- centr absc))))
 
264
                (setf res43
 
265
                        (+ res43 (* fval (f2cl-lib:fref w43b (k) ((1 12))))))
 
266
                (f2cl-lib:fset (f2cl-lib:fref savfun (ipx) ((1 21))) fval)
 
267
               label40))
 
268
            (setf result (* res43 hlgth))
 
269
            (setf abserr
 
270
                    (coerce (abs (* (- res43 res21) hlgth)) 'double-float))
 
271
            (go label65)
 
272
           label45
 
273
            (setf res87 (* (f2cl-lib:fref w87b (23) ((1 23))) fcentr))
 
274
            (setf neval 87)
 
275
            (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 
276
                          ((> k 21) nil)
 
277
              (tagbody
 
278
                (setf res87
 
279
                        (+ res87
 
280
                           (* (f2cl-lib:fref savfun (k) ((1 21)))
 
281
                              (f2cl-lib:fref w87a (k) ((1 21))))))
 
282
               label50))
 
283
            (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 
284
                          ((> k 22) nil)
 
285
              (tagbody
 
286
                (setf absc (* hlgth (f2cl-lib:fref x4 (k) ((1 22)))))
 
287
                (setf res87
 
288
                        (+ res87
 
289
                           (* (f2cl-lib:fref w87b (k) ((1 23)))
 
290
                              (+ (funcall f (+ absc centr))
 
291
                                 (funcall f (- centr absc))))))
 
292
               label60))
 
293
            (setf result (* res87 hlgth))
 
294
            (setf abserr
 
295
                    (coerce (abs (* (- res87 res43) hlgth)) 'double-float))
 
296
           label65
 
297
            (if (and (/= resasc 0.0) (/= abserr 0.0))
 
298
                (setf abserr
 
299
                        (* resasc
 
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))
 
305
           label70))
 
306
       label80
 
307
        (xermsg "SLATEC" "DQNG" "ABNORMAL RETURN" ier 0)
 
308
       label999
 
309
        (go end_label)
 
310
       end_label
 
311
        (return (values nil nil nil nil nil result abserr neval ier))))))
327
312