← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/hypgeo.lisp

  • Committer: Bazaar Package Importer
  • Author(s): Barry deFreese
  • Date: 2006-07-06 17:04:52 UTC
  • mfrom: (1.1.4 upstream)
  • Revision ID: james.westby@ubuntu.com-20060706170452-j9ypoqc1kjfnz221
Tags: 5.9.3-1ubuntu1
* Re-sync with Debian
* Comment out backward-delete-char-untabify in maxima.el (Closes Malone #5273)
* debian/control: build-dep automake -> automake1.9 (Closes BTS: #374663)

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/hypgeo.lisp
2
2
 
3
3
 
4
4
;;    ** (c) Copyright 1976, 1983 Massachusetts Institute of Technology **
5
 
(in-package "MAXIMA")
 
5
(in-package :maxima)
6
6
 
7
7
;;These are the main routines for finding the Laplace Transform
8
8
;; of special functions   --- written by Yannis Avgoustis
113
113
(defun gminc (a b)
114
114
  (list '($gammaincomplete) a b))
115
115
 
 
116
;; Lommel's little s[u,v](z) function.
116
117
(defun littleslommel
117
118
    (m n z)
118
119
  (list '(mqapply)(list '($%s array) m n) z))
469
470
        ((coeffpp)(a zerp)))
470
471
      nil))
471
472
 
472
 
;; Recognize gammagreek(w1,w2), th incomplete gamma function,
 
473
;; Recognize gammagreek(w1,w2), the incomplete gamma function,
473
474
;; integrate(t^(v-1)*exp(-t),t,0,x)
474
475
(defun onegammagreek
475
476
    (exp)
503
504
        ((coeffpp)(a zerp)))
504
505
      nil))
505
506
 
 
507
 
 
508
;; Recognize Lommel s[v1,v2](w) function.
506
509
(defun ones
507
510
    (exp)
508
511
  (m2 exp
513
516
        ((coeffpp)(a zerp)))
514
517
      nil))
515
518
 
 
519
;; Lommel S[v1,v2](w) function
516
520
(defun oneslommel
517
521
    (exp)
518
522
  (m2 exp
927
931
;; Laplace transform of c*t^v*exp(-p*t+e*f).  L contains the pattern
928
932
;; for c*t^v.
929
933
(defun lt-exp (l e f)
930
 
  (prog(c v)
931
 
     (setq c (cdras 'c l) v (cdras 'v l))
932
 
     (cond ((t^2 f)
933
 
            (setq e (inv (mul -8 e)) v (add v 1))
934
 
            (return (f24p146test c v e))))
935
 
     (cond ((sqroott f)
936
 
            (setq e (mul* e e (inv 4)) v (add v 1))
937
 
            (return (f35p147test c v e))))
938
 
     (cond ((t^-1 f)
939
 
            (setq e (mul -4 e) v (add v 1))
940
 
            (return (f29p146test v e))))
941
 
     (return 'other-lt-exponential-to-follow)))
 
934
  (let ((c (cdras 'c l))
 
935
        (v (cdras 'v l)))
 
936
    (cond ((t^2 f)
 
937
           (setq e (inv (mul -8 e)) v (add v 1))
 
938
           (f24p146test c v e))
 
939
          ((sqroott f)
 
940
           (setq e (mul* e e (inv 4)) v (add v 1))
 
941
           (f35p147test c v e))
 
942
          ((t^-1 f)
 
943
           (setq e (mul -4 e) v (add v 1))
 
944
           (f29p146test v e))
 
945
          ((e^-t f)
 
946
           (f36p147 e))
 
947
          ((e^t f)
 
948
           (f37p147 e))
 
949
          (t 'other-lt-exponential-to-follow))))
942
950
 
943
951
(defun t^2(exp)(m2 exp '((mexpt)(t varp) 2) nil))
944
952
 
946
954
 
947
955
(defun t^-1(exp)(m2 exp '((mexpt)(t varp) -1) nil))
948
956
 
 
957
(defun e^-t (exp)
 
958
  (m2 exp
 
959
      '((mexpt) $%e ((mtimes) -1 (t varp)))
 
960
      nil))
 
961
 
 
962
(defun e^t (exp)
 
963
  (m2 exp
 
964
      '((mexpt) $%e (t varp))
 
965
      nil))
 
966
 
949
967
;; Check if conditions for f24p146 hold
950
968
(defun f24p146test (c v a)
951
969
  (cond ((and (eq (asksign a) '$positive)
1021
1039
 
1022
1040
;; Table of Integral Transforms
1023
1041
;;
 
1042
;; p. 147, formula 36:
 
1043
;;
 
1044
;; exp(-a*exp(-t))
 
1045
;;   -> a^(-p)*gamma(p,a)
 
1046
(defun f36p147 (a)
 
1047
  (let ((-a (mul -1 a)))
 
1048
    (mul* (power -a (mul -1 *par*))
 
1049
          `(($gammagreek) ,*par* ,-a))))
 
1050
 
 
1051
;; Table of Integral Transforms
 
1052
;;
 
1053
;; p. 147, formula 36:
 
1054
;;
 
1055
;; exp(-a*exp(t))
 
1056
;;   -> a^(-p)*gammaincomplete(-p,a)
 
1057
(defun f37p147 (a)
 
1058
  (let ((-a (mul -1 a)))
 
1059
    (mul* (power -a *par*)
 
1060
          `(($gammaincomplete) ,(mul -1 *par*) ,-a))))
 
1061
 
 
1062
 
 
1063
;; Table of Integral Transforms
 
1064
;;
1024
1065
;; p. 146, formula 29:
1025
1066
;;
1026
1067
;; t^(v-1)*exp(-a/t/4)
1679
1720
                         ((eq flg 'd)
1680
1721
                          (dtw i1 a1))
1681
1722
                         ((eq flg 'kbateman)
1682
 
                          (kbatemantw a1))
 
1723
                          (kbatemantw i1 a1))
1683
1724
                         ((eq flg 'gammaincomplete)
1684
1725
                          (gammaincompletetw a1 i1))
1685
1726
                         ((eq flg 'kti)
1732
1773
 
1733
1774
(defun sendexec(r a)(distrexecinit ($expand (mul (init r) a)))) 
1734
1775
 
1735
 
(defun whittest
1736
 
    (r a i1 i2)
1737
 
  (cond ((whittindtest i1 i2) 'formula-for-confl-needed)
1738
 
        (t (distrexecinit ($expand (mul (init r)
1739
 
                                        (wtm a i1 i2)))))))
 
1776
;; Test for Whittaker W function.  Simplify this if possible, or
 
1777
;; convert to Whittaker M function.
 
1778
;;
 
1779
;; We have r * %w[i1,i2](a)
 
1780
(defun whittest (r a i1 i2)
 
1781
  (cond ((f16p217test r a i1 i2))
 
1782
        (t
 
1783
         ;; Convert to M function and try again.
 
1784
         (distrexecinit ($expand (mul (init r)
 
1785
                                      (wtm a i1 i2)))))))
1740
1786
 
 
1787
;; Formula 16, p. 217
 
1788
;;
 
1789
;; t^(v-1)*%w[k,u](a*t)
 
1790
;;   -> gamma(u+v+1/2)*gamma(v-u+1/2)*a^(u+1/2)/
 
1791
;;          (gamma(v-k+1)*(p+a/2)^(u+v+1/2)
 
1792
;;        *2f1(u+v+1/2,u-k+1/2;v-k+1;(p-a/2)/(p+a/2))
 
1793
;;
 
1794
;; For Re(v +/- mu) > -1/2
 
1795
(defun f16p217test (r a i1 i2)
 
1796
  ;; We have r*%w[i1,i2](a)
 
1797
  (let ((l (c*t^v r)))
 
1798
    ;; Make sure r is of the form c*t^v
 
1799
    (when l
 
1800
      (let* ((v (add (cdras 'v l) 1))
 
1801
             (c (cdras 'c l)))
 
1802
        ;; Check that v + i2 + 1/2 > 0 and v - i2 + 1/2 > 0.
 
1803
        (when (and (eq (asksign (add (add v i2) 1//2)) '$positive)
 
1804
                   (eq (asksign (add (sub v i2) 1//2)) '$positive))
 
1805
          ;; Ok, we satisfy the conditions.  Now extract the arg.
 
1806
          (let ((l (m2 a
 
1807
                       '((mplus)
 
1808
                         ((coeffpt) (f hasvar) (a freevar))
 
1809
                         ((coeffpp) (c zerp)))
 
1810
                       nil)))
 
1811
            (when l
 
1812
              (let ((a (cdras 'a l)))
 
1813
                ;; We're ready now to compute the transform.
 
1814
                (mul* c
 
1815
                      (power a (add i2 1//2))
 
1816
                      (gm (add (add v i2) 1//2))
 
1817
                      (gm (add (sub v i2) 1//2))
 
1818
                      (inv (mul* (gm (add (sub v i1) 1))
 
1819
                                 (power (add *par* (div a 2))
 
1820
                                        (add (add i2 v) 1//2))))
 
1821
                      (hgfsimp-exec (list (add (add i2 v 1//2))
 
1822
                                          (add (sub i2 i1) 1//2))
 
1823
                                    (list (add (sub v i1) 1))
 
1824
                                    (div (sub *par* (div a 2))
 
1825
                                         (add *par* (div a 2)))))))))))))
 
1826
  
1741
1827
(defun whittindtest (i1 i2)
1742
1828
  (or (maxima-integerp (add i2 i2))
1743
1829
      (neginp (sub (sub (1//2) i2) i1))
1839
1925
           (parcyl (mul* (power 2 inv2) x) -1)))
1840
1926
   (1//2)))
1841
1927
 
1842
 
(defun eitgammaincomplete(x)(mul* -1 (gminc 0 (mul -1 x))))
 
1928
;; The exponential integral Ei can be written in terms of the
 
1929
;; incomplete gamma function.
 
1930
;;
 
1931
;; See Table of Integral Transforms, p. 386:
 
1932
;;
 
1933
;; -Ei(-x) = E1(x) = integrate(exp(-t)/t,t,x,inf)
 
1934
;;
 
1935
;;         = gammaincomplete(0,x)
 
1936
;;
 
1937
(defun eitgammaincomplete (x)
 
1938
  (mul* -1 (gminc 0 (mul -1 x))))
1843
1939
 
 
1940
;; Express Lommel S function in terms of J and Y.
 
1941
;; Luke gives
 
1942
;;
 
1943
;; S[u,v](z) = s[u,v](z) + {2^(u-1)*gamma((u-v+1)/2)*gamma((u+v+1)/2)}
 
1944
;;                 * {sin[(u-v)*%pi/2]*bessel_j(v,z)
 
1945
;;                     - cos[(u-v)*%pi/2]*bessel_y(v,z)
1844
1946
(defun slommeltjandy
1845
1947
    (m n z)
1846
1948
  ((lambda(arg)
1915
2017
;;
1916
2018
;; See Table of Integral Transforms, p.386:
1917
2019
;;
1918
 
;; D[v](z) = 2^(v/2+1/2)*z^(-1/2)*W[v/2+1/4,1/4](z^2/2)
 
2020
;; D[v](z) = 2^(v/2+1/4)*z^(-1/2)*W[v/2+1/4,1/4](z^2/2)
1919
2021
;;
1920
2022
(defun dtw (i a)
1921
 
  (mul* (power 2 (add (div i 2)(inv 4)))
 
2023
  (mul* (power 2 (add (div i 2) (inv 4)))
1922
2024
        (power a (inv -2))
1923
2025
        (wwhit (mul* a a (1//2))
1924
 
               (add (div i 2)(inv 4))
 
2026
               (add (div i 2) (inv 4))
1925
2027
               (inv 4))))
1926
2028
 
1927
2029
;; Bateman's function as a Whittaker W function
1930
2032
;;
1931
2033
;; k[2*v](z) = 1/gamma(v+1)*W[v,1/2](2*z)
1932
2034
;;
1933
 
;; (But this function seems to have v = 1/2.)
1934
 
(defun kbatemantw (a)
1935
 
  ((lambda(ind)
1936
 
     (div (wwhit (add a a) ind (1//2))
1937
 
          (gm (add ind 1))))
1938
 
   (div 1 2)))
 
2035
(defun kbatemantw (v a)
 
2036
  (div (wwhit (add a a) (div v 2) (1//2))
 
2037
       (gm (add (div v 2) 1))))
1939
2038
 
1940
2039
;; Bessel K in terms of Bessel I.
1941
2040
;;
2032
2131
         (div (numjory v sort z 'y)
2033
2132
              (sin% (mul v '$%pi))))))
2034
2133
 
2035
 
;;; LT<foo> functions are various experts on Laplace transforms of the function <foo>.
 
2134
;;; LT<foo> functions are various experts on Laplace transforms of the
 
2135
;;; function <foo>.  The expression being transformed is
 
2136
;;; r*<foo>(args).  The first arg of each expert is r, The remaining
 
2137
;;; args are the arg(s) and/or parameters.
2036
2138
 
2037
 
;;expert on l.t. expressions containing one bessel function of the first kind
 
2139
;; Expert on Laplace transform expressions containing one bessel
 
2140
;; function of the first kind
2038
2141
 
2039
2142
(defun lt1j (rest arg index)
2040
2143
  (lt-ltp 'onej rest arg index))
2042
2145
(defun lt1y (rest arg index)
2043
2146
  (lt-ltp 'oney rest arg index))
2044
2147
 
 
2148
;; Transform of a product of Bessel J functions.  The argument of each
 
2149
;; must be the same, but the orders may be different.
2045
2150
(defun lt2j (rest arg1 arg2 index1 index2)
2046
2151
  (cond ((not (equal arg1 arg2))
2047
2152
         'product-of-bessel-with-different-args)
2050
2155
                   arg1
2051
2156
                   (list 'list index1 index2)))))
2052
2157
 
 
2158
;; Transform of a square of a Bessel J function
2053
2159
(defun lt1j^2
2054
2160
    (rest arg index)
2055
2161
  (lt-ltp 'twoj rest arg (list 'list index index)))
2056
2162
 
 
2163
;; Transform of incomplete Gamma function
2057
2164
(defun lt1gammagreek
2058
2165
    (rest arg1 arg2)
2059
2166
  (lt-ltp 'gammagreek rest arg2 arg1))
2060
2167
 
2061
 
(defun lt1m(r a i1 i2)(lt-ltp 'onem r a (list i1 i2)))
 
2168
;; Laplace transform of r*%m[i1,i2](a)
 
2169
(defun lt1m (r a i1 i2)
 
2170
  (lt-ltp 'onem r a (list i1 i2)))
2062
2171
 
2063
2172
(defun lt1p(r a i1 i2)(lt-ltp 'hyp-onep r a (list i1 i2)))
2064
2173
 
2169
2278
                       (div (add* m n 3) 2))
2170
2279
                 (mul* (inv -4) z z))))
2171
2280
 
2172
 
(defun lt-ltp
2173
 
    (flg rest arg index)
2174
 
  (prog(index1 index2 argl const l l1)
2175
 
     (cond ((or (zerp index)
2176
 
                (eq flg 'onerf)
2177
 
                (eq flg 'onekelliptic)
2178
 
                (eq flg 'onee)
2179
 
                (eq flg 'onepjac)
2180
 
                (eq flg 'd)
2181
 
                (eq flg 's)
2182
 
                (eq flg 'hs)
2183
 
                (eq flg 'ls)
2184
 
                (eq flg 'onem)
2185
 
                (eq flg 'oneq)
2186
 
                (eq flg 'gammagreek)
2187
 
                (eq flg 'asin)
2188
 
                (eq flg 'atan))
2189
 
            (go labl)))
2190
 
     (cond ((or (eq flg 'hyp-onep)(eq flg 'onelog))
 
2281
;; FLG = special function we're transforming
 
2282
;; REST = other stuff
 
2283
;; ARG = arg of special function
 
2284
;; INDEX = index of special function.
 
2285
;;
 
2286
;; So we're transforming REST*FLG(INDEX, ARG).
 
2287
(defun lt-ltp (flg rest arg index)
 
2288
  (prog (index1 index2 argl const l l1)
 
2289
     (when (or (zerp index)
 
2290
               (eq flg 'onerf)
 
2291
               (eq flg 'onekelliptic)
 
2292
               (eq flg 'onee)
 
2293
               (eq flg 'onepjac)
 
2294
               (eq flg 'd)
 
2295
               (eq flg 's)
 
2296
               (eq flg 'hs)
 
2297
               (eq flg 'ls)
 
2298
               (eq flg 'onem)
 
2299
               (eq flg 'oneq)
 
2300
               (eq flg 'gammagreek)
 
2301
               (eq flg 'asin)
 
2302
               (eq flg 'atan))
 
2303
            (go labl))
 
2304
     (cond ((or (eq flg 'hyp-onep)
 
2305
                (eq flg 'onelog))
 
2306
            ;; Go to labl1 if we've got %p or log.
2191
2307
            (go labl1)))
2192
 
     (cond ((not (consp index)) (go lab)))
2193
 
     (cond ((not (eq (car index) 'list))(go lab)))
 
2308
     ;; Skip this if we have exactly one index or if INDEX doesn't
 
2309
     ;; look like '(LIST i1 i2)
 
2310
     (cond ((not (consp index))
 
2311
            (go lab)))
 
2312
     (cond ((not (eq (car index) 'list))
 
2313
            (go lab)))
 
2314
     ;; If the first index is exactly 0, set INDEX1 to it and go to
 
2315
     ;; LA.
2194
2316
     (cond ((zerp (setq index1 (cadr index)))(go la)))
 
2317
     ;; If we're here, the index is of the form '(LIST i1 i2), which
 
2318
     ;; means this is the product of two Bessel functions.
 
2319
     ;;
 
2320
     ;; Ok.  I think this is wrong, in general.  I think we get here
 
2321
     ;; if we're doing the produce to two Bessel functions.  This code
 
2322
     ;; seems to be applying the property that bessel_j(-n,x) =
 
2323
     ;; (-1)^n*bessel_j(n,x).  But this only works if n is an integer.
 
2324
     #+nil
2195
2325
     (cond ((eq (checksigntm (simplifya (inv (setq index1
2196
2326
                                                   (cadr
2197
2327
                                                    index)))
2198
2328
                                        nil))
2199
2329
                '$negative)
 
2330
            ;; FIXME: What is this supposed to do?  We take the
 
2331
            ;; reciprocal of the first index and see if it's negative.
 
2332
            ;; If so, we change the sign of index1 and divide REST by
 
2333
            ;; the index.  What is this for?
2200
2334
            (setq index1
2201
2335
                  (mul -1 index1)
2202
2336
                  rest
2203
2337
                  (mul* (power -1 index1) rest))))
2204
2338
     la
2205
 
     (cond ((zerp (setq index2 (caddr index)))(go la2)))
 
2339
     ;; If the second index is zero, skip over this.
 
2340
     (cond ((zerp (setq index2 (caddr index)))
 
2341
            (go la2)))
 
2342
     ;; Wrong too.  See comment above about bessel_j(-n,x).
 
2343
     #+nil
2206
2344
     (cond ((eq (checksigntm (simplifya (inv (setq index2
2207
2345
                                                   (caddr
2208
2346
                                                    index)))
2209
2347
                                        nil))
2210
2348
                '$negative)
 
2349
            ;; FIXME: This does the same for index2 as for index1
 
2350
            ;; above.  Why?
2211
2351
            (setq index2
2212
2352
                  (mul -1 index2)
2213
2353
                  rest
2214
2354
                  (mul* (power -1 index2) rest))))
2215
2355
     la2
 
2356
     ;; Put the 2 indices in a list and go on.
2216
2357
     (setq index (list index1 index2))
2217
2358
     (go labl)
2218
2359
     lab
 
2360
     ;; We're here if we have one index, and it's not one of the
 
2361
     ;; special cases.
 
2362
     ;;
 
2363
     ;; FIXME:  Find out what functions trigger this.
 
2364
 
 
2365
     ;; I think this is the one bessel function case with a negative
 
2366
     ;; index.  At least here we check that the index is an integer.
 
2367
     ;; We can either leave it here or take it out.  It seems
 
2368
     ;; bessel_j(-n,x) is converted by the Bessel simplifiers before
 
2369
     ;; we get here.
2219
2370
     (cond ((and (eq (checksigntm (simplifya (inv index)
2220
2371
                                             nil))
2221
2372
                     '$negative)
2222
2373
                 (maxima-integerp index))
 
2374
            ;; FIXME: Same as the 2 index thing above, but we need an
 
2375
            ;; (numeric) integer too.  Why?
2223
2376
            (setq index (mul -1 index))
2224
2377
            (setq rest (mul (power -1 index) rest))))
2225
2378
     labl
 
2379
     ;; Handle index = 0 or one of the special functions erf,
 
2380
     ;; kelliptic, E, Jacobi, %d, %s, hstruve, lstruve, %m, Q,
 
2381
     ;; incomplete gamma, asin, atan.
2226
2382
     (setq argl (f+c arg))
2227
 
     (setq const (cdras 'c argl) arg (cdras 'f argl))
 
2383
     (setq const (cdras 'c argl)
 
2384
           arg (cdras 'f argl))
 
2385
     ;; See if the arg is f + c, and replace arg with f.
2228
2386
     (cond ((null const)(go labl1)))
 
2387
     ;; This handles the case of when the const term is actually there.
2229
2388
     (cond ((not (eq (checksigntm (simplifya (power const
2230
2389
                                                    2)
2231
2390
                                             nil))
2232
2391
                     '$zero))
 
2392
            ;; I guess prop4 handles the case when square of the
 
2393
            ;; constant term is not zero.  Too bad I (rtoy) don't know
 
2394
            ;; what prop4 is.
 
2395
            ;;
 
2396
            ;; FIXME:  Implement prop4.
 
2397
            (format t "const = ~A~%" const)
 
2398
            (format t "f = ~A~%" arg)
2233
2399
            (return 'prop4-to-be-applied)))
2234
2400
     labl1
2235
 
     (cond ((eq flg 'oney)(return (lty rest arg index))))
 
2401
     ;; No const term, if we're here.
 
2402
     (cond ((eq flg 'oney)
 
2403
            ;; Handle bessel_y here.  We're done.
 
2404
            (return (lty rest arg index))))
 
2405
     ;; Try to express the function in terms of hypergeometric
 
2406
     ;; functions that we can handle.
2236
2407
     (cond ((setq l
2237
2408
                  (d*x^m*%e^a*x ($factor (mul* rest
2238
2409
                                               (car (setq
2241
2412
                                                      flg
2242
2413
                                                      index
2243
2414
                                                      arg)))))))
 
2415
            ;; Convert the special function to a hypgergeometric
 
2416
            ;; function.  L1 is the special function converted to the
 
2417
            ;; hypergeometric function.  d*x^m*%e^a*x looks for that
 
2418
            ;; factor in the expanded form.
2244
2419
            (return (%$etest l l1))))
 
2420
     ;; We currently don't know how to handle this yet.
2245
2421
     (return 'other-ca-later)))
2246
2422
 
2247
2423
(defun lty
2312
2488
    (asin (asintf arg))
2313
2489
    (atan (atantf arg))))
2314
2490
 
 
2491
;; Create a hypergeometric form that we recognize.  The function is
 
2492
;; %f[n,m](p; q; arg).  We represent this as a list of the form
 
2493
;; (fpq (<length p> <length q>) <p> <q> <arg>)
2315
2494
(defun ref-fpq (p q arg)
2316
2495
  (list 'fpq (list (length p) (length q))
2317
2496
        p q arg))
2496
2675
                 (list (add 1 n) (add 1 m) (add* 1 n m))
2497
2676
                 (mul -1 (power arg 2)))))
2498
2677
 
2499
 
;; Match d*x^m*%e^(a*x)
 
2678
;; Match d*x^m*%e^(a*x).  If we match, Q is the e^(a*x) part, A is a,
 
2679
;; M is M, and D is d.
2500
2680
(defun d*x^m*%e^a*x
2501
2681
    (exp)
2502
2682
  (m2 exp
2523
2703
     (cond ((eq (car ans) 'dionimo)
2524
2704
            (return (dionarghyp-y l index (cadr ans)))))
2525
2705
     (return 'fail-in-execfy)))
2526
 
;;executive  for recognizing the sort of argument
2527
2706
 
2528
 
(defun execargmatch
2529
 
    (arg)
 
2707
;; Executive for recognizing the sort of argument to the
 
2708
;; hypergeometric function.  We look to see if the arg is of the form
 
2709
;; a*x^m + c.  Return a list of 'dionimo (what does that mean?) and
 
2710
;; the match.
 
2711
(defun execargmatch (arg)
2530
2712
  (prog(l1)
2531
2713
     (cond ((setq l1 (a*x^m+c ($factor arg)))
2532
2714
            (return (list 'dionimo l1))))
2534
2716
            (return (list 'dionimo l1))))
2535
2717
     (return 'other-case-args-to-follow)))
2536
2718
 
 
2719
;; We have hypergeometric function whose arg looks like a*x^m+c.  L1
 
2720
;; matches the d*x^m... part, L2 is the hypergeometric function and
 
2721
;; arg is the match for a*x^m+c.
2537
2722
(defun dionarghyp
2538
2723
    (l1 l2 arg)
2539
2724
  (prog(a m c)
2547
2732
            (return (f19cond a m l1 l2))))
2548
2733
     (return 'prop4-and-aother-cases-to-folow)))
2549
2734
 
2550
 
 
 
2735
 
 
2736
;; Match f(x)+c
2551
2737
(defun f+c
2552
2738
    (exp)
2553
2739
  (m2 exp
2554
2740
      '((mplus)((coeffpt)(f hasvar))((coeffpp)(c freevar)))
2555
2741
      nil))
2556
2742
 
 
2743
;; Match a*x^m+c.
2557
2744
(defun a*x^m+c
2558
2745
    (exp)
2559
2746
  (m2 exp
2598
2785
           d (cdras 'd l1)
2599
2786
           l1 (caddr l2)
2600
2787
           l2 (cadddr l2))
 
2788
     ;; At this point, we have the function d*x^s*%f[p,q](l1, l2, (a*t)^m).
 
2789
     ;; Check to see if Formula 19, p 220 applies.
2601
2790
     (cond ((and (not (eq (checksigntm (sub (add* p
2602
2791
                                                  m
2603
2792
                                                  -1)