← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/cpoly.lisp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Sauthier
  • Date: 2009-07-13 15:38:41 UTC
  • mfrom: (3.1.3 squeeze)
  • Revision ID: james.westby@ubuntu.com-20090713153841-gtux06oun30kuuo7
Tags: 5.17.1-1ubuntu1
https://launchpad.net/bugs/296643
https://launchpad.net/bugs/242243
* Merge from debian unstable, remaining changes (LP: #296643, LP: #242243):
   - debian/maxima-doc.doc-base.{tips, plotting}:
    + Use .shtml instead of .html to fix lintian errors.
   - debian/maxima-emacs.emacsen-install:
    + Install symlinks for source files rather than copying them.  This
      makes find-function work.
    + Install symlink for *.lisp so that we don't need to add
      /usr/share/emacs/site-lisp/maxima to load-path.
  - debian/maxima-emacs.emacsen-startup:
    + Remove use of /usr/share/emacs/site-lisp/maxima, since this
      causes load-path shadows and is not needed anymore.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/cpoly.lisp
84
84
                      (return (simplify (list '(mplus)
85
85
                                              (simplify (list '(mtimes) '$%i ci))
86
86
                                              cr))))
87
 
                    0d0 0d0))
 
87
                    0.0 0.0))
88
88
                  (t (return (list '(mlist simp)))))))
89
89
     (setq degree (cadr expr) nn (1+ degree))
90
 
     (setq *pr-sl* (make-array nn :initial-element 0d0))
91
 
     (setq *pi-sl* (make-array nn :initial-element 0d0))
 
90
     (setq *pr-sl* (make-array nn :initial-element 0.0))
 
91
     (setq *pi-sl* (make-array nn :initial-element 0.0))
92
92
     (or (catch 'notpoly
93
93
           (errset (do ((expr (cdr expr) (cddr expr)) (l) (%i (cadr res)))
94
94
                       ((null expr))
104
104
                            (setq complex t))))))
105
105
         ;;this should catch expressions like sin(x)-x
106
106
         (cpoly-err expr1))
107
 
     (setq *shr-sl* (make-array nn :initial-element 0d0))
108
 
     (setq *shi-sl* (make-array nn :initial-element 0d0))
109
 
     (setq *qpr-sl* (make-array nn :initial-element 0d0))
110
 
     (setq *hr-sl*  (make-array degree :initial-element 0d0))
111
 
     (setq *qhr-sl* (make-array degree :initial-element 0d0))
 
107
     (setq *shr-sl* (make-array nn :initial-element 0.0))
 
108
     (setq *shi-sl* (make-array nn :initial-element 0.0))
 
109
     (setq *qpr-sl* (make-array nn :initial-element 0.0))
 
110
     (setq *hr-sl*  (make-array degree :initial-element 0.0))
 
111
     (setq *qhr-sl* (make-array degree :initial-element 0.0))
112
112
     (when complex
113
 
       (setq *qpi-sl* (make-array nn :initial-element 0d0))
114
 
       (setq *hi-sl*  (make-array degree :initial-element 0d0))
115
 
       (setq *qhi-sl* (make-array degree :initial-element 0d0)))
 
113
       (setq *qpi-sl* (make-array nn :initial-element 0.0))
 
114
       (setq *hi-sl*  (make-array degree :initial-element 0.0))
 
115
       (setq *qhi-sl* (make-array degree :initial-element 0.0)))
116
116
     (setq nn degree)
117
117
     (if complex
118
118
         (setq res (errset (cpoly-sl degree)))
124
124
     (setq res nil)
125
125
     (cond ((not (zerop nn))
126
126
            (mtell "~%Only ~S out of ~S roots found " (- degree nn) degree)
127
 
            (setq expr 0d0)
 
127
            (setq expr 0.0)
128
128
            (do ((i 0 (1+ i)))
129
129
                ((> i nn))
130
130
              (setq expr
142
142
                                     (float (car den))
143
143
                                     (float (cadr den)))
144
144
                          (simplify (list '(mplus) (simplify (list '(mtimes) '$%i ci)) cr)))
145
 
                        0d0 0d0)
 
145
                        0.0 0.0)
146
146
                  res (cons expr res))))
147
147
     (do ((i degree (1- i)))
148
148
         ((= i nn))
172
172
(defun cpoly-sl (degree) 
173
173
  ((lambda (logbas infin smalno are mre xx yy cosr sinr cr ci sr si tr
174
174
            ti zr zi bnd n polysc polysc1 conv) 
175
 
     (setq mre (* 2d0 (sqrt 2d0) are)
 
175
     (setq mre (* 2.0 (sqrt 2.0) are)
176
176
           yy (- xx))
177
177
     (do ((i degree (1- i)))
178
178
         ((not (and (zerop (aref *pr-sl* i)) (zerop (aref *pi-sl* i))))
182
182
     (do ((i 0 (1+ i)))
183
183
         ((> i nn))
184
184
       (setf (aref *shr-sl* i) (cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
185
 
     (scale-sl)
 
185
     (if (> nn 0) (scale-sl))
186
186
     (do ()
187
187
         ((> 2 nn)
188
188
          (cdivid-sl (- (aref *pr-sl* 1)) (- (aref *pi-sl* 1))
224
224
       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))
225
225
       (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) j)))
226
226
     nn)
227
 
   (log 2d0)
228
 
   most-positive-double-float
229
 
   least-positive-double-float
230
 
   double-float-epsilon
231
 
   0d0 0.70710677d0 0d0 -0.069756474d0 0.99756405d0
232
 
   0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0 0 0 nil))
 
227
   (log 2.0)
 
228
   most-positive-flonum
 
229
   least-positive-flonum
 
230
   flonum-epsilon
 
231
   0.0
 
232
   ;; sqrt(0.5)
 
233
   0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207863675069231154561485124624180279253686063220607485499679157066113329637527963778999752505763910302857350547799858029851372672984310073642587093204445993047761646152421543571607254198813018139976257039948436267
 
234
   0.0
 
235
   ;; cos(94deg)
 
236
   -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804
 
237
   ;; sin(94deg)
 
238
   0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315
 
239
   0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0 0 nil))
233
240
 
234
241
(defun noshft-sl (l1) 
235
242
  (do ((i 0 (1+ i))
241
248
  (do ((jj 1 (1+ jj)))
242
249
      ((> jj l1))
243
250
    (cond ((> (cmod-sl (aref *hr-sl* n) (aref *hi-sl* n))
244
 
              (* 10d0 are (cmod-sl (aref *pr-sl* n) (aref *pi-sl* n))))
 
251
              (* 10.0 are (cmod-sl (aref *pr-sl* n) (aref *pi-sl* n))))
245
252
           (cdivid-sl (- (aref *pr-sl* nn)) (- (aref *pi-sl* nn)) (aref *hr-sl* n) (aref *hi-sl* n))
246
253
           (setq tr cr ti ci)
247
254
           (do ((j n (1- j)) (t1) (t2))
255
262
                 ((> 1 j))
256
263
               (setf (aref *hr-sl* j) (aref *hr-sl* (1- j)))
257
264
               (setf (aref *hi-sl* j) (aref *hi-sl* (1- j))))
258
 
             (setf (aref *hr-sl* 0) 0d0)
259
 
             (setf (aref *hi-sl* 0) 0d0))))) 
 
265
             (setf (aref *hr-sl* 0) 0.0)
 
266
             (setf (aref *hi-sl* 0) 0.0))))) 
260
267
 
261
268
(defun fxshft-sl (l2) 
262
269
  ((lambda (test pasd otr oti svsr svsi bool pvr pvi) 
291
298
                    ((setq pasd nil))))))
292
299
     (or conv (vrshft-sl 10))
293
300
     nil)
294
 
   t nil 0d0 0d0 0d0 0d0 nil 0d0 0d0)) 
 
301
   t nil 0.0 0.0 0.0 0.0 nil 0.0 0.0)) 
295
302
 
296
303
(defun vrshft-sl (l3) 
297
304
  (setq conv nil sr zr si zi)
301
308
      ((> i l3))
302
309
    (polyev-sl)
303
310
    (setq mp (cmod-sl pvr pvi) ms (cmod-sl sr si))
304
 
    (cond ((> (* 20d0 (errev-sl ms mp)) mp)
 
311
    (cond ((> (* 20.0 (errev-sl ms mp)) mp)
305
312
           (setq conv t zr sr zi si)
306
313
           (return t)))
307
314
    (cond ((= i 1) (setq omp mp))
308
315
          ((or bool1 (> omp mp) (not (< relstp 0.05)))
309
 
           (if (> (* 0.1d0 mp) omp)
 
316
           (if (> (* 0.1 mp) omp)
310
317
               (return t)
311
318
               (setq omp mp)))
312
319
          (t (setq tp relstp bool1 t)
332
339
       (hvr (setf (aref *qhr-sl* 0) (aref *hr-sl* 0)))
333
340
       (hvi (setf (aref *qhi-sl* 0) (aref *hi-sl* 0))))
334
341
      ((> i n)
335
 
       (setq bool (not (> (cmod-sl hvr hvi) (* 10d0 are (cmod-sl (aref *hr-sl* n) (aref *hi-sl* n))))))
 
342
       (setq bool (not (> (cmod-sl hvr hvi) (* 10.0 are (cmod-sl (aref *hr-sl* n) (aref *hi-sl* n))))))
336
343
       (cond ((not bool) (cdivid-sl (- pvr) (- pvi) hvr hvi) (setq tr cr ti ci))
337
 
             (t (setq tr 0d0 ti 0d0)))
 
344
             (t (setq tr 0.0 ti 0.0)))
338
345
       nil)
339
346
    (setq $t (- (+ (aref *hr-sl* i) (* hvr sr)) (* hvi si)))
340
347
    (setf (aref *qhi-sl* i) (setq hvi (+ (aref *hi-sl* i) (* hvr si) (* hvi sr))))
345
352
                  ((> j n))
346
353
                (setf (aref *hr-sl* j) (aref *qhr-sl* (1- j)))
347
354
                (setf (aref *hi-sl* j) (aref *qhi-sl* (1- j))))
348
 
              (setf (aref *hr-sl* 0) 0d0)
349
 
              (setf (aref *hi-sl* 0) 0d0))
 
355
              (setf (aref *hr-sl* 0) 0.0)
 
356
              (setf (aref *hi-sl* 0) 0.0))
350
357
        (t (do ((j 1. (1+ j)) (t1) (t2))
351
358
               ((> j n))
352
359
             (setq t1 (aref *qhr-sl* (1- j)) t2 (aref *qhi-sl* (1- j)))
381
388
            (cond ((> x xm) (setq x xm)))))
382
389
     (do ((f))
383
390
         (nil)
384
 
       (setq xm (* 0.1d0 x) f (aref *shr-sl* 0))
 
391
       (setq xm (* 0.1 x) f (aref *shr-sl* 0))
385
392
       (do ((i 1 (1+ i))) ((> i nn)) (setq f (+ (aref *shr-sl* i) (* f xm))))
386
 
       (cond ((not (< 0d0 f)) (return t)))
 
393
       (cond ((not (< 0.0 f)) (return t)))
387
394
       (setq x xm))
388
395
     (do ((dx x) (df) (f))
389
 
         ((> 5d-3 (abs (/ dx x))) x)
 
396
         ((> 5e-3 (abs (/ dx x))) x)
390
397
       (setq f (aref *shr-sl* 0) df f)
391
398
       (do ((i 1 (1+ i)))
392
399
           ((> i n))
393
400
         (setq f (+ (* f x) (aref *shr-sl* i)) df (+ (* df x) f)))
394
401
       (setq f (+ (* f x) (aref *shr-sl* nn)) dx (/ f df) x (- x dx))))
395
402
   (exp (/ (- (log (aref *shr-sl* nn)) (log (aref *shr-sl* 0))) (float nn)))
396
 
   0d0)) 
 
403
   0.0)) 
397
404
 
398
405
(defun scale-sl ()
399
 
  (do ((i 0 (1+ i)) (j 0) (x 0d0) (dx 0d0))
 
406
  (do ((i 0 (1+ i)) (j 0) (x 0.0) (dx 0.0))
400
407
      ((> i nn)
401
408
       (setq x (/ x (float (- (1+ nn) j)))
402
409
             dx (/ (- (log (aref *shr-sl* nn)) (log (aref *shr-sl* 0))) (float nn))
403
 
             polysc1 (floor (+ 0.5d0 (/ dx logbas)))
404
 
             x (+ x (* (float (* polysc1 nn)) logbas 0.5d0))
405
 
             polysc (floor (+ 0.5d0 (/ x logbas)))))
 
410
             polysc1 (floor (+ 0.5 (/ dx logbas)))
 
411
             x (+ x (* (float (* polysc1 nn)) logbas 0.5))
 
412
             polysc (floor (+ 0.5 (/ x logbas)))))
406
413
    (cond ((zerop (aref *shr-sl* i)) (setq j (1+ j)))
407
414
          (t (setq x (+ x (log (aref *shr-sl* i)))))))
408
415
  (do ((i nn (1- i)) (j (- polysc) (+ j polysc1)))
426
433
                  cr (/ br bi) 
427
434
                  br (- ai (* ar r1)) 
428
435
                  ci (/ br bi)))))
429
 
   0d0)
 
436
   0.0)
430
437
  nil) 
431
438
 
432
439
(defun cmod-sl (ar ai) 
434
441
        ai (abs ai))
435
442
  (cond ((> ai ar) (setq ar (/ ar ai)) (* ai (sqrt (1+ (* ar ar)))))
436
443
        ((> ar ai) (setq ai (/ ai ar)) (* ar (sqrt (1+ (* ai ai)))))
437
 
        ((* 1.41421357d0 ar)))) 
 
444
        ((* 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872534 ; sqrt(2.0)
 
445
          ar))))
438
446
 
439
 
;;;this is the algorithm for doing real polynomials.
440
 
;;;it is algorithm 493 from acm toms vol 1 p 178 (1975) by jenkins.
441
 
;;;note that array indexing starts from 0.
442
 
;;;the names of the arrays have been changed to be the same as for cpoly.
443
 
;;;the correspondence is:  p - pr-sl, qp - qpr-sl, k - hr-sl, qk - qhr-sl,
444
 
;;;svk - shr-sl, temp - shi-sl.  the roots are put in pr-sl and pi-sl.
445
 
;;;the variable si appears not to be used here
 
447
;;; This is the algorithm for doing real polynomials.
 
448
;;; It is algorithm 493 from acm toms vol 1 p 178 (1975) by jenkins.
 
449
;;; Note that array indexing starts from 0.
 
450
;;; The names of the arrays have been changed to be the same as for cpoly.
 
451
;;; The correspondence is:  p - pr-sl, qp - qpr-sl, k - hr-sl, qk - qhr-sl,
 
452
;;; svk - shr-sl, temp - shi-sl.  the roots are put in pr-sl and pi-sl.
 
453
;;; The variable si appears not to be used here
446
454
 
447
455
(declare-top (special sr u v a b c d a1 a3 a7 e f g h szr szi lzr lzi are
448
456
                      mre n nn nz ui vi s $polyfactor)) 
454
462
     (do ((i degree (1- i))) ((not (zerop (aref *pr-sl* i))) (setq nn i n (1- i))))
455
463
     (setq degree nn)
456
464
     (do ((i 0 (1+ i))) ((> i nn)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
457
 
     (scale-sl)
 
465
     (if (> nn 0) (scale-sl))
458
466
     (do nil
459
467
         ((< nn 3)
460
468
          (cond ((= nn 2)
462
470
                 (cond ((and $polyfactor (not (zerop szi)))
463
471
                        (setf (aref *pr-sl* 2) (/ (aref *pr-sl* 2) (aref *pr-sl* 0)))
464
472
                        (setf (aref *pr-sl* 1) (/ (aref *pr-sl* 1) (aref *pr-sl* 0)))
465
 
                        (setf (aref *pi-sl* 2) 1d0))
 
473
                        (setf (aref *pi-sl* 2) 1.0))
466
474
                       (t (setf (aref *pr-sl* 2) szr)
467
475
                          (setf (aref *pi-sl* 2) szi)
468
476
                          (setf (aref *pr-sl* 1) lzr)
482
490
         (cond (zerok (do ((j n (1- j)))
483
491
                          ((< j 1))
484
492
                        (setf (aref *hr-sl* j) (aref *hr-sl* (1- j))))
485
 
                      (setf (aref *hr-sl* 0) 0d0)
 
493
                      (setf (aref *hr-sl* 0) 0.0)
486
494
                      (setq zerok (zerop (aref *hr-sl* n))))
487
495
               (t (setq t1 (- (/ aa cc)))
488
496
                  (do ((j n (1- j)))
490
498
                    (setf (aref *hr-sl* j) (+ (* t1 (aref *hr-sl* (1- j))) (aref *pr-sl* j))))
491
499
                  (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
492
500
                  (setq zerok (not (> (abs (aref *hr-sl* n))
493
 
                                      (* (abs bb) are 10d0)))))))
 
501
                                      (* (abs bb) are 10.0)))))))
494
502
       (do ((i 0 (1+ i))) ((> i n)) (setf (aref *shi-sl* i) (aref *hr-sl* i)))
495
503
       (do ((cnt 1 (1+ cnt)))
496
504
           ((> cnt 20.) (setq conv1 nil))
498
506
                      (- (* cosr xx) (* sinr yy))
499
507
                    (setq yy (+ (* sinr xx) (* cosr yy)))) 
500
508
               sr (* bnd xx) 
501
 
               u (* -2d0 sr) 
 
509
               u (* -2.0 sr) 
502
510
               v bnd)
503
511
         (fxshfr-sl (* 20 cnt))
504
512
         (cond ((> nz 0)
510
518
                       (cond ((and $polyfactor (not (zerop szi)))
511
519
                              (setf (aref *pr-sl* nn) v)
512
520
                              (setf (aref *pr-sl* n) u)
513
 
                              (setf (aref *pi-sl* nn) 1d0)))))
 
521
                              (setf (aref *pi-sl* nn) 1.0)))))
514
522
                (setq nn (- nn nz) n (1- nn))
515
523
                (do ((i 0 (1+ i))) ((> i nn)) (setf (aref *pr-sl* i) (aref *qpr-sl* i)))
516
524
                (return nil)))
531
539
     (do ((i 0 (1+ i)) (j (- polysc (* polysc1 degree)) (+ j polysc1)))
532
540
         ((> i nn))
533
541
       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))))
534
 
   (log 2d0)
535
 
   most-positive-double-float
536
 
   least-positive-double-float
537
 
   double-float-epsilon
538
 
   0d0 0.70710677d0 0d0 -0.069756474d0 0.99756405d0
539
 
   0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0 0 0 0 0 t))  
 
542
   (log 2.0)
 
543
   most-positive-flonum
 
544
   least-positive-flonum
 
545
   flonum-epsilon
 
546
   0.0
 
547
   ;; sqrt(0.5)
 
548
   0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207863675069231154561485124624180279253686063220607485499679157066113329637527963778999752505763910302857350547799858029851372672984310073642587093204445993047761646152421543571607254198813018139976257039948436267
 
549
   0.0
 
550
   ;; cos(94deg)
 
551
   -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804
 
552
   ;; sin(94deg)
 
553
   0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315
 
554
   0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0 0 0 0 t))  
540
555
 
541
556
(defun fxshfr-sl (l2) 
542
557
  ((lambda (my-type a b c d e f g h a1 a3 a7)
545
560
     (quadsd-sl)
546
561
     (calcsc-sl)
547
562
     (do ((j 1 (1+ j))
548
 
          (betav 0.25d0)
549
 
          (betas 0.25d0)
 
563
          (betav 0.25)
 
564
          (betas 0.25)
550
565
          (oss sr)
551
566
          (ovv v)
552
567
          (tvv) (tss) (ss) (vv) (tv) (ts) (ots) (otv)
556
571
       (calcsc-sl)
557
572
       (newest-sl)
558
573
       (setq vv vi
559
 
             ss 0d0)
 
574
             ss 0.0)
560
575
       (or (zerop (aref *hr-sl* n))
561
576
           (setq ss (- (/ (aref *pr-sl* nn) (aref *hr-sl* n)))))
562
 
       (setq tv 1d0 ts 1d0)
 
577
       (setq tv 1.0 ts 1.0)
563
578
       (cond ((not (or (= j 1) (= my-type 3)))
564
579
              (or (zerop vv) (setq tv (abs (/ (- vv ovv) vv))))
565
580
              (or (zerop ss) (setq ts (abs (/ (- ss oss) ss))))
566
 
              (setq tvv 1d0)
 
581
              (setq tvv 1.0)
567
582
              (and (< tv otv) (setq tvv (* tv otv)))
568
 
              (setq tss 1d0)
 
583
              (setq tss 1.0)
569
584
              (and (< ts ots) (setq tss (* ts ots)))
570
585
              (setq vpass (< tvv betav) spass (< tss betas))
571
586
              (cond ((or spass vpass)
580
595
                            (cond (bool (quadit-sl)
581
596
                                        (and (> nz 0) (return t))
582
597
                                        (setq vtry t
583
 
                                              betav (* 0.25d0 betav))
 
598
                                              betav (* 0.25 betav))
584
599
                                        (cond ((or stry (not spass))
585
600
                                               (setq l50 t))
586
601
                                              (t (do ((i 0 (1+ i)))
590
605
                            (cond ((not l50)
591
606
                                   (setq iflag (realit-sl))
592
607
                                   (and (> nz 0) (return t))
593
 
                                   (setq stry t betas (* 0.25d0 betas))
 
608
                                   (setq stry t betas (* 0.25 betas))
594
609
                                   (cond ((zerop iflag) (setq l50 t))
595
610
                                         (t (setq ui (- (+ s s)) 
596
611
                                                  vi (* s s))))))
608
623
             oss ss
609
624
             otv tv
610
625
             ots ts)))
611
 
   0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0))
 
626
   0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0))
612
627
 
613
628
(defun quadit-sl nil 
614
629
  (setq nz 0 u ui v vi)
615
630
  (do ((tried) (j 0) (ee) (zm) (t1) (mp) (relstp) (omp))
616
631
      (nil)
617
 
    (quad-sl 1d0 u v)
618
 
    (and (> (abs (- (abs szr) (abs lzr))) (* 1d-2 (abs lzr)))
 
632
    (quad-sl 1.0 u v)
 
633
    (and (> (abs (- (abs szr) (abs lzr))) (* 1e-2 (abs lzr)))
619
634
         (return nil))
620
635
    (quadsd-sl)
621
636
    (setq mp (+ (abs (- a (* szr b))) (abs (* szi b))) 
622
637
          zm (sqrt (abs v)) 
623
 
          ee (* 2d0 (abs (aref *qpr-sl* 0)))
 
638
          ee (* 2.0 (abs (aref *qpr-sl* 0)))
624
639
          t1 (- (* szr b)))
625
640
    (do ((i 1 (1+ n)))
626
641
        ((> i n)) (setq ee (+ (* ee zm) (abs (aref *qpr-sl* i)))))
627
642
    (setq ee (+ (* ee zm) (abs (+ a t1))) 
628
 
          ee (- (* (+ (* 5d0 mre) (* 4d0 are)) ee)
629
 
                 (* (+ (* 5d0 mre) (* 2d0 are))
 
643
          ee (- (* (+ (* 5.0 mre) (* 4.0 are)) ee)
 
644
                 (* (+ (* 5.0 mre) (* 2.0 are))
630
645
                     (+ (abs (+ a t1)) (* (abs b) zm)))
631
 
                 (* -2d0 are (abs t1))))
632
 
    (cond ((not (> mp (* 2d1 ee))) (setq nz 2)
 
646
                 (* -2.0 are (abs t1))))
 
647
    (cond ((not (> mp (* 2e1 ee))) (setq nz 2)
633
648
           (return nil)))
634
649
    (setq j (1+ j))
635
650
    (and (> j 20) (return nil))
636
 
    (cond ((not (or (< j 2) (> relstp 1d-2) (< mp omp) tried))
 
651
    (cond ((not (or (< j 2) (> relstp 1e-2) (< mp omp) tried))
637
652
           (and (< relstp are) (setq relstp are))
638
653
           (setq relstp (sqrt relstp) 
639
654
                 u (- u (* u relstp)) 
665
680
          ee (* (/ mre (+ are mre)) (abs (aref *qpr-sl* 0))))
666
681
    (do ((i 1 (1+ i)))
667
682
        ((> i nn)) (setq ee (+ (* ee ms) (abs (aref *qpr-sl* i)))))
668
 
    (cond ((not (> mp (* 2d1 (- (* (+ are mre) ee) (* mre mp)))))
669
 
           (setq nz 1 szr s szi 0d0)
 
683
    (cond ((not (> mp (* 2e1 (- (* (+ are mre) ee) (* mre mp)))))
 
684
           (setq nz 1 szr s szi 0.0)
670
685
           (return 0)))
671
686
    (setq j (1+ j))
672
687
    (and (> j 10) (return 0))
673
688
    (cond ((not (or (< j 2)
674
 
                    (> (abs t1) (* 1d-3 (abs (- s t1))))
 
689
                    (> (abs t1) (* 1e-3 (abs (- s t1))))
675
690
                    (not (> mp omp))))
676
691
           (return 1)))
677
692
    (setq omp mp kv (aref *hr-sl* 0))
680
695
        ((> i n))
681
696
      (setq kv (+ (* kv s) (aref *hr-sl* i)))
682
697
      (setf (aref *qhr-sl* i) kv))
683
 
    (cond ((> (abs kv) (* (abs (aref *hr-sl* n)) 1d1 are))
 
698
    (cond ((> (abs kv) (* (abs (aref *hr-sl* n)) 1e1 are))
684
699
           (setq t1 (- (/ pv kv)))
685
700
           (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
686
701
           (do ((i 1 (1+ i)))
687
702
               ((> i n))
688
703
             (setf (aref *hr-sl* i)
689
704
                    (+ (* t1 (aref *qhr-sl* (1- i))) (aref *qpr-sl* i)))))
690
 
          (t (setf (aref *hr-sl* 0) 0d0)
 
705
          (t (setf (aref *hr-sl* 0) 0.0)
691
706
             (do ((i 1 (1+ i)))
692
707
                 ((> i n)) (setf (aref *hr-sl* i) (aref *qhr-sl* (1- i))))))
693
708
    (setq kv (aref *hr-sl* 0))
694
709
    (do ((i 1 (1+ i)))
695
710
        ((> i n)) (setq kv (+ (* kv s) (aref *hr-sl* i))))
696
 
    (setq t1 0d0)
697
 
    (and (> (abs kv) (* (abs (aref *hr-sl* n)) 10d0 are))
 
711
    (setq t1 0.0)
 
712
    (and (> (abs kv) (* (abs (aref *hr-sl* n)) 10.0 are))
698
713
         (setq t1 (- (/ pv kv))))
699
714
    (setq s (+ s t1))))
700
715
 
710
725
    (setq c0 (- (aref *hr-sl* i) (* u c) (* v d)))
711
726
    (setf (aref *qhr-sl* i) c0)
712
727
    (setq d c c c0))
713
 
  (cond ((not (or (> (abs c) (* (abs (aref *hr-sl* n)) 1d2 are))
714
 
                  (> (abs d) (* (abs (aref *hr-sl* (1- n))) 1d2 are))))
 
728
  (cond ((not (or (> (abs c) (* (abs (aref *hr-sl* n)) 1e2 are))
 
729
                  (> (abs d) (* (abs (aref *hr-sl* (1- n))) 1e2 are))))
715
730
         (setq my-type 3))
716
731
        ((not (< (abs d) (abs c)))
717
732
         (setq my-type 2 
735
750
(defun nextk-sl ()
736
751
  (declare (special my-type))
737
752
  (cond ((= my-type 3)
738
 
         (setf (aref *hr-sl* 0) 0d0)
739
 
         (setf (aref *hr-sl* 1) 0d0)
 
753
         (setf (aref *hr-sl* 0) 0.0)
 
754
         (setf (aref *hr-sl* 1) 0.0)
740
755
         (do ((i 2 (1+ i)))
741
756
             ((> i n)) (setf (aref *hr-sl* i) (aref *qhr-sl* (- i 2)))))
742
 
        ((> (abs a1) (* (abs (if (= my-type 1) b a)) 1d1 are))
 
757
        ((> (abs a1) (* (abs (if (= my-type 1) b a)) 1e1 are))
743
758
         (setq a7 (/ a7 a1) a3 (/ a3 a1))
744
759
         (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
745
760
         (setf (aref *hr-sl* 1) (- (aref *qpr-sl* 1) (* a7 (aref *qpr-sl* 0))))
749
764
                  (+ (* a3 (aref *qhr-sl* (- i 2)))
750
765
                      (- (* a7 (aref *qpr-sl* (1- i))))
751
766
                      (aref *qpr-sl* i)))))
752
 
        (t (setf (aref *hr-sl* 0) 0d0)
 
767
        (t (setf (aref *hr-sl* 0) 0.0)
753
768
           (setf (aref *hr-sl* 1) (- (* a7 (aref *qpr-sl* 0))))
754
769
           (do ((i 2 (1+ i)))
755
770
               ((> i n))
761
776
(defun newest-sl ()
762
777
  (declare (special my-type))
763
778
  ((lambda (a4 a5 b1 b2 c1 c2 c3 c4) 
764
 
     (cond ((= my-type 3) (setq ui 0d0 vi 0d0))
 
779
     (cond ((= my-type 3) (setq ui 0.0 vi 0.0))
765
780
           (t (if (= my-type 2)
766
781
                  (setq a4 (+ (* (+ a g) f) h) 
767
782
                        a5 (+ (* (+ f u) c) (* v d)))
775
790
                    c4 (- c1 c2 c3) 
776
791
                    c1 (+ a5 (* b1 a4) (- c4)))
777
792
              (if (zerop c1)
778
 
                  (setq ui 0d0 vi 0d0)
 
793
                  (setq ui 0.0 vi 0.0)
779
794
                  (setq ui (- u (/ (+ (* u (+ c3 c2))
780
795
                                      (* v (+ (* b1 a1) (* b2 a7))))
781
796
                                   c1))
782
797
                        vi (* v (1+ (/ c4 c1)))))))
783
798
     nil)
784
 
   0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0)) 
 
799
   0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0)) 
785
800
 
786
801
(defun quadsd-sl ()
787
802
  (setq b (aref *pr-sl* 0))
797
812
          a c0)))
798
813
 
799
814
(defun quad-sl (a0 b1 c0) 
800
 
  (setq szr 0d0 szi 0d0 lzr 0d0 lzi 0d0)
801
 
  ((lambda (b0 d0 e) 
 
815
  (setq szr 0.0 szi 0.0 lzr 0.0 lzi 0.0)
 
816
  ((lambda (b0 l0 e) 
802
817
     (cond ((zerop a0) (or (zerop b1) (setq szr (- (/ c0 b1)))))
803
818
           ((zerop c0) (setq lzr (- (/ b1 a0))))
804
 
           (t (setq b0 (/ b1 2d0))
 
819
           (t (setq b0 (/ b1 2.0))
805
820
              (cond ((< (abs b0) (abs c0))
806
821
                     (setq e a0)
807
 
                     (and (< c0 0.0d0) (setq e (- a0)))
 
822
                     (and (< c0 0.0) (setq e (- a0)))
808
823
                     (setq e (- (* b0 (/ b0 (abs c0))) e) 
809
 
                           d0 (* (sqrt (abs e)) (sqrt (abs c0)))))
810
 
                    (t (setq e (- 1d0 (* (/ a0 b0) (/ c0 b0))) 
811
 
                             d0 (* (sqrt (abs e)) (abs b0)))))
812
 
              (cond ((< e 0.0d0)
 
824
                           l0 (* (sqrt (abs e)) (sqrt (abs c0)))))
 
825
                    (t (setq e (- 1.0 (* (/ a0 b0) (/ c0 b0))) 
 
826
                             l0 (* (sqrt (abs e)) (abs b0)))))
 
827
              (cond ((< e 0.0)
813
828
                     (setq szr (- (/ b0 a0)) 
814
829
                           lzr szr 
815
 
                           szi (abs (/ d0 a0)) 
 
830
                           szi (abs (/ l0 a0)) 
816
831
                           lzi (- szi)))
817
 
                    (t (or (< b0 0d0) (setq d0 (- d0)))
818
 
                       (setq lzr (/ (- d0 b0) a0))
 
832
                    (t (or (< b0 0.0) (setq l0 (- l0)))
 
833
                       (setq lzr (/ (- l0 b0) a0))
819
834
                       (or (zerop lzr) (setq szr (/ c0 lzr a0)))))))
820
835
     nil)
821
 
   0d0 0d0 0d0)) 
 
836
   0.0 0.0 0.0)) 
822
837
 
823
838
(declare-top (unspecial logbas infin smalno are mre cr ci sr si tr ti zr zi
824
839
                        n nn bool conv pvr pvi polysc polysc1 sr u v