← Back to branch summary

~ubuntu-branches/ubuntu/karmic/maxima/karmic

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2004-11-13 18:39:14 UTC
  • mto: (2.1.2 hoary) (3.2.1 sid) (1.1.5 upstream)
  • mto: This revision was merged to the branch mainline in revision 3.
  • Revision ID: james.westby@ubuntu.com-20041113183914-ttig0evwuatnqosl
Tags: upstream-5.9.1
Import upstream version 5.9.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/numerical/slatec/dqk15i.lisp
 
1
;;; Compiled by f2cl version 2.0 beta 2002-05-06
 
2
;;; 
 
3
;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
 
4
;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
 
5
;;;           (:array-slicing t) (:declare-common nil)
 
6
;;;           (:float-format double-float))
 
7
 
 
8
(in-package "SLATEC")
 
9
 
 
10
 
 
11
(let ((wg (make-array 8 :element-type 'double-float))
 
12
      (xgk (make-array 8 :element-type 'double-float))
 
13
      (wgk (make-array 8 :element-type 'double-float)))
 
14
  (declare (type (simple-array double-float (8)) wgk xgk wg))
 
15
  (f2cl-lib:fset (f2cl-lib:fref wg (1) ((1 8))) 0.0)
 
16
  (f2cl-lib:fset (f2cl-lib:fref wg (2) ((1 8))) 0.1294849661688697)
 
17
  (f2cl-lib:fset (f2cl-lib:fref wg (3) ((1 8))) 0.0)
 
18
  (f2cl-lib:fset (f2cl-lib:fref wg (4) ((1 8))) 0.27970539148927664)
 
19
  (f2cl-lib:fset (f2cl-lib:fref wg (5) ((1 8))) 0.0)
 
20
  (f2cl-lib:fset (f2cl-lib:fref wg (6) ((1 8))) 0.3818300505051189)
 
21
  (f2cl-lib:fset (f2cl-lib:fref wg (7) ((1 8))) 0.0)
 
22
  (f2cl-lib:fset (f2cl-lib:fref wg (8) ((1 8))) 0.4179591836734694)
 
23
  (f2cl-lib:fset (f2cl-lib:fref xgk (1) ((1 8))) 0.9914553711208126)
 
24
  (f2cl-lib:fset (f2cl-lib:fref xgk (2) ((1 8))) 0.9491079123427585)
 
25
  (f2cl-lib:fset (f2cl-lib:fref xgk (3) ((1 8))) 0.8648644233597691)
 
26
  (f2cl-lib:fset (f2cl-lib:fref xgk (4) ((1 8))) 0.7415311855993945)
 
27
  (f2cl-lib:fset (f2cl-lib:fref xgk (5) ((1 8))) 0.5860872354676911)
 
28
  (f2cl-lib:fset (f2cl-lib:fref xgk (6) ((1 8))) 0.4058451513773972)
 
29
  (f2cl-lib:fset (f2cl-lib:fref xgk (7) ((1 8))) 0.20778495500789848)
 
30
  (f2cl-lib:fset (f2cl-lib:fref xgk (8) ((1 8))) 0.0)
 
31
  (f2cl-lib:fset (f2cl-lib:fref wgk (1) ((1 8))) 0.022935322010529224)
 
32
  (f2cl-lib:fset (f2cl-lib:fref wgk (2) ((1 8))) 0.06309209262997856)
 
33
  (f2cl-lib:fset (f2cl-lib:fref wgk (3) ((1 8))) 0.10479001032225019)
 
34
  (f2cl-lib:fset (f2cl-lib:fref wgk (4) ((1 8))) 0.14065325971552592)
 
35
  (f2cl-lib:fset (f2cl-lib:fref wgk (5) ((1 8))) 0.1690047266392679)
 
36
  (f2cl-lib:fset (f2cl-lib:fref wgk (6) ((1 8))) 0.19035057806478542)
 
37
  (f2cl-lib:fset (f2cl-lib:fref wgk (7) ((1 8))) 0.20443294007529889)
 
38
  (f2cl-lib:fset (f2cl-lib:fref wgk (8) ((1 8))) 0.20948214108472782)
 
39
  (defun dqk15i (f boun inf a b result abserr resabs resasc)
 
40
    (declare (type f2cl-lib:integer4 inf)
 
41
     (type double-float resasc resabs abserr result b a boun)
 
42
     (type (function (double-float) (values double-float &rest t)) f))
 
43
    (prog ((fv1 (make-array 7 :element-type 'double-float))
 
44
           (fv2 (make-array 7 :element-type 'double-float)) (j 0) (absc 0.0)
 
45
           (absc1 0.0) (absc2 0.0) (centr 0.0) (dinf 0.0) (epmach 0.0) (fc 0.0)
 
46
           (fsum 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (resg 0.0) (resk 0.0)
 
47
           (reskh 0.0) (tabsc1 0.0) (tabsc2 0.0) (uflow 0.0) (abs$ 0.0f0))
 
48
      (declare (type single-float abs$)
 
49
       (type (simple-array double-float (7)) fv2 fv1)
 
50
       (type double-float uflow tabsc2 tabsc1 reskh resk resg hlgth fval2 fval1
 
51
        fsum fc epmach dinf centr absc2 absc1 absc)
 
52
       (type f2cl-lib:integer4 j))
 
53
      (setf epmach (f2cl-lib:d1mach 4))
 
54
      (setf uflow (f2cl-lib:d1mach 1))
 
55
      (setf dinf
 
56
              (coerce
 
57
               (the f2cl-lib:integer4
 
58
                    (min (the f2cl-lib:integer4 1)
 
59
                         (the f2cl-lib:integer4 inf)))
 
60
               'double-float))
 
61
      (setf centr (* 0.5 (+ a b)))
 
62
      (setf hlgth (* 0.5 (- b a)))
 
63
      (setf tabsc1 (+ boun (/ (* dinf (- 1.0 centr)) centr)))
 
64
      (setf fval1
 
65
              (multiple-value-bind
 
66
                  (ret-val var-0)
 
67
                  (funcall f tabsc1)
 
68
                (declare (ignore))
 
69
                (when var-0 (setf tabsc1 var-0))
 
70
                ret-val))
 
71
      (if (= inf 2) (setf fval1 (+ fval1 (funcall f (- tabsc1)))))
 
72
      (setf fc (/ (/ fval1 centr) centr))
 
73
      (setf resg (* (f2cl-lib:fref wg (8) ((1 8))) fc))
 
74
      (setf resk (* (f2cl-lib:fref wgk (8) ((1 8))) fc))
 
75
      (setf resabs (coerce (abs resk) 'double-float))
 
76
      (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 
77
                    ((> j 7) nil)
 
78
        (tagbody
 
79
          (setf absc (* hlgth (f2cl-lib:fref xgk (j) ((1 8)))))
 
80
          (setf absc1 (- centr absc))
 
81
          (setf absc2 (+ centr absc))
 
82
          (setf tabsc1 (+ boun (/ (* dinf (- 1.0 absc1)) absc1)))
 
83
          (setf tabsc2 (+ boun (/ (* dinf (- 1.0 absc2)) absc2)))
 
84
          (setf fval1
 
85
                  (multiple-value-bind
 
86
                      (ret-val var-0)
 
87
                      (funcall f tabsc1)
 
88
                    (declare (ignore))
 
89
                    (when var-0 (setf tabsc1 var-0))
 
90
                    ret-val))
 
91
          (setf fval2
 
92
                  (multiple-value-bind
 
93
                      (ret-val var-0)
 
94
                      (funcall f tabsc2)
 
95
                    (declare (ignore))
 
96
                    (when var-0 (setf tabsc2 var-0))
 
97
                    ret-val))
 
98
          (if (= inf 2) (setf fval1 (+ fval1 (funcall f (- tabsc1)))))
 
99
          (if (= inf 2) (setf fval2 (+ fval2 (funcall f (- tabsc2)))))
 
100
          (setf fval1 (/ (/ fval1 absc1) absc1))
 
101
          (setf fval2 (/ (/ fval2 absc2) absc2))
 
102
          (f2cl-lib:fset (f2cl-lib:fref fv1 (j) ((1 7))) fval1)
 
103
          (f2cl-lib:fset (f2cl-lib:fref fv2 (j) ((1 7))) fval2)
 
104
          (setf fsum (+ fval1 fval2))
 
105
          (setf resg (+ resg (* (f2cl-lib:fref wg (j) ((1 8))) fsum)))
 
106
          (setf resk (+ resk (* (f2cl-lib:fref wgk (j) ((1 8))) fsum)))
 
107
          (setf resabs
 
108
                  (+ resabs
 
109
                     (* (f2cl-lib:fref wgk (j) ((1 8)))
 
110
                        (+ (abs fval1) (abs fval2)))))
 
111
         label10))
 
112
      (setf reskh (* resk 0.5))
 
113
      (setf resasc (* (f2cl-lib:fref wgk (8) ((1 8))) (abs (- fc reskh))))
 
114
      (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 
115
                    ((> j 7) nil)
 
116
        (tagbody
 
117
          (setf resasc
 
118
                  (+ resasc
 
119
                     (* (f2cl-lib:fref wgk (j) ((1 8)))
 
120
                        (+ (abs (- (f2cl-lib:fref fv1 (j) ((1 7))) reskh))
 
121
                           (abs (- (f2cl-lib:fref fv2 (j) ((1 7))) reskh))))))
 
122
         label20))
 
123
      (setf result (* resk hlgth))
 
124
      (setf resasc (* resasc hlgth))
 
125
      (setf resabs (* resabs hlgth))
 
126
      (setf abserr (coerce (abs (* (- resk resg) hlgth)) 'double-float))
 
127
      (if (and (/= resasc 0.0) (/= abserr 0.0))
 
128
          (setf abserr
 
129
                  (* resasc (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
 
130
      (if (> resabs (/ uflow (* 50.0 epmach)))
 
131
          (setf abserr (max (* epmach 50.0 resabs) abserr)))
 
132
      (go end_label)
 
133
     end_label
 
134
      (return (values nil nil nil nil nil result abserr resabs resasc)))))
 
135