?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Remove the obsolete DEFMTRFUN-EXTERNAL macro
[maxima.git] / share / hompack / lisp / fixpnf.lisp
blob826c11895b151a4a2fa3240e3988a99d29442a86
1 ;;; Compiled by f2cl version:
2 ;;; ("f2cl1.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
3 ;;; "f2cl2.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
4 ;;; "f2cl3.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
5 ;;; "f2cl4.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
6 ;;; "f2cl5.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
7 ;;; "f2cl6.l,v 1d5cbacbb977 2008/08/24 00:56:27 rtoy $"
8 ;;; "macros.l,v 1409c1352feb 2013/03/24 20:44:50 toy $")
9
10 ;;; Using Lisp CMU Common Lisp snapshot-2020-04 (21D Unicode)
11 ;;;
12 ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
13 ;;; (:coerce-assigns :as-needed) (:array-type ':array)
14 ;;; (:array-slicing t) (:declare-common nil)
15 ;;; (:float-format double-float))
16
17 (in-package "HOMPACK")
18
19
20 (let* ((limitd 1000) (cursw (coerce 10.0f0 'double-float)))
21 (declare (type (f2cl-lib:integer4 1000 1000) limitd)
22 (type (double-float) cursw)
23 (ignorable limitd cursw))
24 (let ((s 0.0)
25 (relerr 0.0)
26 (hold 0.0)
27 (h 0.0)
28 (curtol 0.0)
29 (abserr 0.0)
30 (np1 0)
31 (nfec 0)
32 (nc 0)
33 (limit 0)
34 (jw 0)
35 (iter 0)
36 (iflagc 0)
37 (start nil)
38 (polsys nil)
39 (crash nil))
40 (declare (type (double-float) s relerr hold h curtol abserr)
41 (type (f2cl-lib:integer4) np1 nfec nc limit jw iter iflagc)
42 (type f2cl-lib:logical start polsys crash))
43 (labels ((multi-entry-fixpnf
44 (%name% n y iflag arcre arcae ansre ansae trace$ a nfe arclen
45 yp yold ypold qr alpha tz pivot w wp z0 z1 sspar par ipar)
46 (declare (type (array f2cl-lib:integer4 (*)) ipar)
47 (type (array double-float (*)) par)
48 (type (array double-float (*)) sspar)
49 (type (array f2cl-lib:integer4 (*)) pivot)
50 (type (double-float) arclen ansae ansre arcae arcre)
51 (type (array double-float (*)) z1 z0 wp w tz alpha qr
52 ypold yold yp a y)
53 (type (f2cl-lib:integer4) nfe trace$ iflag n))
54 (f2cl-lib:with-multi-array-data
55 ((y double-float y-%data% y-%offset%)
56 (a double-float a-%data% a-%offset%)
57 (yp double-float yp-%data% yp-%offset%)
58 (yold double-float yold-%data% yold-%offset%)
59 (ypold double-float ypold-%data% ypold-%offset%)
60 (qr double-float qr-%data% qr-%offset%)
61 (alpha double-float alpha-%data% alpha-%offset%)
62 (tz double-float tz-%data% tz-%offset%)
63 (w double-float w-%data% w-%offset%)
64 (wp double-float wp-%data% wp-%offset%)
65 (z0 double-float z0-%data% z0-%offset%)
66 (z1 double-float z1-%data% z1-%offset%)
67 (pivot f2cl-lib:integer4 pivot-%data% pivot-%offset%)
68 (sspar double-float sspar-%data% sspar-%offset%)
69 (par double-float par-%data% par-%offset%)
70 (ipar f2cl-lib:integer4 ipar-%data% ipar-%offset%))
71 (prog ()
72 (declare)
73 (if (eq %name% 'polynf) (go polynf))
74 (setf polsys f2cl-lib:%false%)
75 (go label11)
76 polynf
77 (setf polsys f2cl-lib:%true%)
78 label11
79 (if (or (<= n 0) (<= ansre 0.0f0) (< ansae 0.0f0))
80 (setf iflag 7))
81 (if (and (>= iflag -2) (<= iflag 0)) (go label20))
82 (if (= iflag 2) (go label120))
83 (if (= iflag 3) (go label90))
84 (setf iflag 7)
85 (go end_label)
86 label20
87 (setf arclen (coerce 0.0f0 'double-float))
88 (if (<= arcre 0.0f0)
89 (setf arcre (* 0.5f0 (f2cl-lib:fsqrt ansre))))
90 (if (<= arcae 0.0f0)
91 (setf arcae (* 0.5f0 (f2cl-lib:fsqrt ansae))))
92 (setf nc n)
93 (setf nfec 0)
94 (setf iflagc iflag)
95 (setf np1 (f2cl-lib:int-add n 1))
96 (setf start f2cl-lib:%true%)
97 (setf crash f2cl-lib:%false%)
98 (setf hold (coerce 1.0f0 'double-float))
99 (setf h (coerce 0.1f0 'double-float))
100 (setf s (coerce 0.0f0 'double-float))
101 (setf (f2cl-lib:fref ypold-%data%
102 (1)
103 ((1 (f2cl-lib:int-add n 1)))
104 ypold-%offset%)
105 (coerce 1.0f0 'double-float))
106 (setf (f2cl-lib:fref yp-%data%
107 (1)
108 ((1 (f2cl-lib:int-add n 1)))
109 yp-%offset%)
110 (coerce 1.0f0 'double-float))
111 (setf (f2cl-lib:fref y-%data%
112 (1)
113 ((1 (f2cl-lib:int-add n 1)))
114 y-%offset%)
115 (coerce 0.0f0 'double-float))
116 (f2cl-lib:fdo (jw 2 (f2cl-lib:int-add jw 1))
117 ((> jw np1) nil)
118 (tagbody
119 (setf (f2cl-lib:fref ypold-%data%
120 (jw)
121 ((1 (f2cl-lib:int-add n 1)))
122 ypold-%offset%)
123 (coerce 0.0f0 'double-float))
124 (setf (f2cl-lib:fref yp-%data%
125 (jw)
126 ((1 (f2cl-lib:int-add n 1)))
127 yp-%offset%)
128 (coerce 0.0f0 'double-float))
129 label40))
130 (if
131 (<= (f2cl-lib:fref sspar-%data% (1) ((1 8)) sspar-%offset%)
132 0.0f0)
133 (setf (f2cl-lib:fref sspar-%data%
134 (1)
135 ((1 8))
136 sspar-%offset%)
137 (coerce 0.5f0 'double-float)))
138 (if
139 (<= (f2cl-lib:fref sspar-%data% (2) ((1 8)) sspar-%offset%)
140 0.0f0)
141 (setf (f2cl-lib:fref sspar-%data%
142 (2)
143 ((1 8))
144 sspar-%offset%)
145 (coerce 0.01f0 'double-float)))
146 (if
147 (<= (f2cl-lib:fref sspar-%data% (3) ((1 8)) sspar-%offset%)
148 0.0f0)
149 (setf (f2cl-lib:fref sspar-%data%
150 (3)
151 ((1 8))
152 sspar-%offset%)
153 (coerce 0.5f0 'double-float)))
154 (if
155 (<= (f2cl-lib:fref sspar-%data% (4) ((1 8)) sspar-%offset%)
156 0.0f0)
157 (setf (f2cl-lib:fref sspar-%data%
158 (4)
159 ((1 8))
160 sspar-%offset%)
161 (* (+ (f2cl-lib:fsqrt (+ n 1.0f0)) 4.0f0)
162 (f2cl-lib:d1mach 4))))
163 (if
164 (<= (f2cl-lib:fref sspar-%data% (5) ((1 8)) sspar-%offset%)
165 0.0f0)
166 (setf (f2cl-lib:fref sspar-%data%
167 (5)
168 ((1 8))
169 sspar-%offset%)
170 (coerce 1.0f0 'double-float)))
171 (if
172 (<= (f2cl-lib:fref sspar-%data% (6) ((1 8)) sspar-%offset%)
173 0.0f0)
174 (setf (f2cl-lib:fref sspar-%data%
175 (6)
176 ((1 8))
177 sspar-%offset%)
178 (coerce 0.1f0 'double-float)))
179 (if
180 (<= (f2cl-lib:fref sspar-%data% (7) ((1 8)) sspar-%offset%)
181 0.0f0)
182 (setf (f2cl-lib:fref sspar-%data%
183 (7)
184 ((1 8))
185 sspar-%offset%)
186 (coerce 3.0f0 'double-float)))
187 (if
188 (<= (f2cl-lib:fref sspar-%data% (8) ((1 8)) sspar-%offset%)
189 0.0f0)
190 (setf (f2cl-lib:fref sspar-%data%
191 (8)
192 ((1 8))
193 sspar-%offset%)
194 (coerce 2.0f0 'double-float)))
195 (cond
196 ((>= iflagc (f2cl-lib:int-sub 1))
197 (f2cl-lib:fdo (jw 2 (f2cl-lib:int-add jw 1))
198 ((> jw np1) nil)
199 (tagbody
200 (setf (f2cl-lib:fref a-%data%
201 ((f2cl-lib:int-sub jw 1))
202 ((1 n))
203 a-%offset%)
204 (f2cl-lib:fref y-%data%
205 (jw)
206 ((1 (f2cl-lib:int-add n 1)))
207 y-%offset%))
208 label60))))
209 label90
210 (setf limit limitd)
211 label120
212 (f2cl-lib:fdo (iter 1 (f2cl-lib:int-add iter 1))
213 ((> iter limit) nil)
214 (tagbody
215 (cond
216 ((< (f2cl-lib:fref y (1) ((1 (f2cl-lib:int-add n 1))))
217 0.0f0)
218 (setf arclen s)
219 (setf iflag 5)
220 (go end_label)))
221 label140
222 (setf curtol (* cursw hold))
223 (setf relerr arcre)
224 (setf abserr arcae)
225 (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
226 ((> jw np1) nil)
227 (tagbody
228 (cond
229 ((>
230 (abs
231 (+
232 (f2cl-lib:fref yp
233 (jw)
234 ((1 (f2cl-lib:int-add n 1))))
235 (-
236 (f2cl-lib:fref ypold
237 (jw)
238 ((1
239 (f2cl-lib:int-add n 1)))))))
240 curtol)
241 (setf relerr ansre)
242 (setf abserr ansae)
243 (go label200)))
244 label160))
245 label200
246 (multiple-value-bind
247 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7
248 var-8 var-9 var-10 var-11 var-12 var-13 var-14
249 var-15 var-16 var-17 var-18 var-19 var-20 var-21
250 var-22 var-23 var-24 var-25)
251 (stepnf nc nfec iflagc start crash hold h relerr
252 abserr s y yp yold ypold a qr alpha tz pivot w wp
253 z0 z1 sspar par ipar)
254 (declare (ignore var-0 var-10 var-11 var-12 var-13
255 var-14 var-15 var-16 var-17 var-18
256 var-19 var-20 var-21 var-22 var-23
257 var-24 var-25))
258 (setf nfec var-1)
259 (setf iflagc var-2)
260 (setf start var-3)
261 (setf crash var-4)
262 (setf hold var-5)
263 (setf h var-6)
264 (setf relerr var-7)
265 (setf abserr var-8)
266 (setf s var-9))
267 (cond
268 ((> trace$ 0)
269 (f2cl-lib:fformat trace
270 ("~%" " STEP" 1 (("~5D")) "~3@T"
271 "NFE =" 1 (("~5D")) "~3@T"
272 "ARC LENGTH =" 1 (("~9,4,0,'*,F"))
273 "~3@T" "LAMBDA =" 1
274 (("~7,4,0,'*,F")) "~5@T"
275 "X vector:" "~%" t
276 ("~1@T" 6 (("~12,4,2,0,'*,,'EE")))
277 "~%")
278 iter
279 nfec
280 s
281 (f2cl-lib:fref y-%data%
282 (1)
283 ((1
284 (f2cl-lib:int-add
285 n
286 1)))
287 y-%offset%)
288 (do ((jw 2 (f2cl-lib:int-add jw 1))
289 (%ret nil))
290 ((> jw np1) (nreverse %ret))
291 (declare (type f2cl-lib:integer4 jw))
292 (push
293 (f2cl-lib:fref y-%data%
294 (jw)
295 ((1
296 (f2cl-lib:int-add
297 n
298 1)))
299 y-%offset%)
300 %ret)))))
301 (setf nfe nfec)
302 (cond
303 ((> iflagc 0)
304 (setf arclen s)
305 (setf iflag iflagc)
306 (go end_label)))
307 (cond
308 (crash
309 (setf iflag 2)
310 (if (< arcre relerr) (setf arcre relerr))
311 (if (< ansre relerr) (setf ansre relerr))
312 (if (< arcae abserr) (setf arcae abserr))
313 (if (< ansae abserr) (setf ansae abserr))
314 (setf limit (f2cl-lib:int-sub limit iter))
315 (go end_label)))
316 (cond
317 ((>=
318 (f2cl-lib:fref y (1) ((1 (f2cl-lib:int-add n 1))))
319 1.0f0)
320 (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
321 ((> jw np1) nil)
322 (tagbody
323 (setf (f2cl-lib:fref z0-%data%
324 (jw)
325 ((1 (f2cl-lib:int-add n 1)))
326 z0-%offset%)
327 (f2cl-lib:fref yold-%data%
328 (jw)
329 ((1
330 (f2cl-lib:int-add n 1)))
331 yold-%offset%))
332 label260))
333 (multiple-value-bind
334 (var-0 var-1 var-2 var-3 var-4 var-5 var-6
335 var-7 var-8 var-9 var-10 var-11 var-12 var-13
336 var-14 var-15 var-16 var-17)
337 (rootnf nc nfec iflagc ansre ansae y yp yold
338 ypold a qr alpha tz pivot w wp par ipar)
339 (declare (ignore var-0 var-3 var-4 var-5 var-6
340 var-7 var-8 var-9 var-10 var-11
341 var-12 var-13 var-14 var-15 var-16
342 var-17))
343 (setf nfec var-1)
344 (setf iflagc var-2))
345 (setf nfe nfec)
346 (setf iflag 1)
347 (if (> iflagc 0) (setf iflag iflagc))
348 (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
349 ((> jw np1) nil)
350 (tagbody
351 (setf (f2cl-lib:fref w-%data%
352 (jw)
353 ((1 (f2cl-lib:int-add n 1)))
354 w-%offset%)
355 (-
356 (f2cl-lib:fref y-%data%
357 (jw)
358 ((1
359 (f2cl-lib:int-add n
360 1)))
361 y-%offset%)
362 (f2cl-lib:fref z0-%data%
363 (jw)
364 ((1
365 (f2cl-lib:int-add n
366 1)))
367 z0-%offset%)))
368 label290))
369 (setf arclen (+ (- s hold) (dnrm2 np1 w 1)))
370 (go end_label)))
371 (cond
372 (polsys
373 (cond
374 ((<
375 (f2cl-lib:fref yp
376 (1)
377 ((1 (f2cl-lib:int-add n 1))))
378 0.0f0)
379 (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
380 ((> jw np1) nil)
381 (tagbody
382 (setf (f2cl-lib:fref yp-%data%
383 (jw)
384 ((1
385 (f2cl-lib:int-add n
386 1)))
387 yp-%offset%)
388 (-
389 (f2cl-lib:fref yp-%data%
390 (jw)
391 ((1
392 (f2cl-lib:int-add n
393 1)))
394 yp-%offset%)))
395 (setf (f2cl-lib:fref ypold-%data%
396 (jw)
397 ((1
398 (f2cl-lib:int-add n
399 1)))
400 ypold-%offset%)
401 (f2cl-lib:fref yp-%data%
402 (jw)
403 ((1
404 (f2cl-lib:int-add n
405 1)))
406 yp-%offset%))
407 label310))
408 (setf start f2cl-lib:%true%)))))
409 label400))
410 (setf iflag 3)
411 (setf arclen s)
412 (go end_label)
413 end_label
414 (return
415 (values nil
416 nil
417 iflag
418 arcre
419 arcae
420 ansre
421 ansae
422 nil
423 nil
424 nfe
425 arclen
426 nil
427 nil
428 nil
429 nil
430 nil
431 nil
432 nil
433 nil
434 nil
435 nil
436 nil
437 nil
438 nil
439 nil))))))
440 (defun fixpnf
441 (n y iflag arcre arcae ansre ansae trace$ a nfe arclen yp yold
442 ypold qr alpha tz pivot w wp z0 z1 sspar par ipar)
443 (multiple-value-bind
444 (v0 v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12 v13 v14 v15 v16 v17
445 v18 v19 v20 v21 v22 v23 v24)
446 (multi-entry-fixpnf 'fixpnf n y iflag arcre arcae ansre ansae
447 trace$ a nfe arclen yp yold ypold qr alpha tz pivot w wp z0 z1
448 sspar par ipar)
449 (values v0
450 v1
451 v2
452 v3
453 v4
454 v5
455 v6
456 v7
457 v8
458 v9
459 v10
460 v11
461 v12
462 v13
463 v14
464 v15
465 v16
466 v17
467 v18
468 v19
469 v20
470 v21
471 v22
472 v23
473 v24)))
474 (defun polynf
475 (n y iflag arcre arcae ansre ansae trace$ a nfe arclen yp yold
476 ypold qr alpha tz pivot w wp z0 z1 sspar par ipar)
477 (multiple-value-bind
478 (v0 v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12 v13 v14 v15 v16 v17
479 v18 v19 v20 v21 v22 v23 v24)
480 (multi-entry-fixpnf 'polynf n y iflag arcre arcae ansre ansae
481 trace$ a nfe arclen yp yold ypold qr alpha tz pivot w wp z0 z1
482 sspar par ipar)
483 (values v0
484 v1
485 v2
486 v3
487 v4
488 v5
489 v6
490 v7
491 v8
492 v9
493 v10
494 v11
495 v12
496 v13
497 v14
498 v15
499 v16
500 v17
501 v18
502 v19
503 v20
504 v21
505 v22
506 v23
507 v24))))))
508
509 (in-package #:cl-user)
510 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
511 (eval-when (:load-toplevel :compile-toplevel :execute)
512 (setf (gethash 'fortran-to-lisp::fixpnf
513 fortran-to-lisp::*f2cl-function-info*)
514 (fortran-to-lisp::make-f2cl-finfo
515 :arg-types '((fortran-to-lisp::integer4) (array double-float (*))
516 (fortran-to-lisp::integer4) (double-float)
517 (double-float) (double-float) (double-float)
518 (fortran-to-lisp::integer4) (array double-float (*))
519 (fortran-to-lisp::integer4) (double-float)
520 (array double-float (*)) (array double-float (*))
521 (array double-float (*)) (array double-float (*))
522 (array double-float (*)) (array double-float (*))
523 (array fortran-to-lisp::integer4 (*))
524 (array double-float (*)) (array double-float (*))
525 (array double-float (*)) (array double-float (*))
526 (array double-float (*)) (array double-float (*))
527 (array fortran-to-lisp::integer4 (*)))
528 :return-values '(nil nil fortran-to-lisp::iflag
529 fortran-to-lisp::arcre fortran-to-lisp::arcae
530 fortran-to-lisp::ansre fortran-to-lisp::ansae nil
531 nil fortran-to-lisp::nfe fortran-to-lisp::arclen
532 nil nil nil nil nil nil nil nil nil nil nil nil nil
533 nil)
534 :calls '(fortran-to-lisp::dnrm2 fortran-to-lisp::rootnf
535 fortran-to-lisp::stepnf fortran-to-lisp::d1mach))))
536
Maxima, a Computer Algebra System