Changeset 0d5870
 Timestamp:
 03/22/12 01:45:44 (3 years ago)
 Branches:
 master
 Children:
 9d3daf
 Parents:
 405df6
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

qdgamma.lisp
rc388f8 r0d5870 480 480 (+ k v 1))))))) 481 481 482 483 ;; For v not an integer: 484 ;; 485 ;; E(v,z) = gamma(1v)*z^(v1)  sum((1)^k*z^k/(kv+1)/k!, k, 0, inf) 486 ;; 487 ;; For v an integer: 488 ;; 489 ;; E(v,z) = (z)^(v1)/(v1)!*(psi(v)log(z)) 490 ;;  sum((1)^k*z^k/(kv+1)/k!, k, 0, inf, k != n1) 491 ;; 482 492 (defun sexpintegrale (v z) 483 493 ;; E(v,z) = gamma(1v)*z^(v1)  sum((1)^k*z^k/(kv+1)/k!, k, 0, inf) … … 485 495 (v ( v)) 486 496 (eps (epsilon z))) 487 (loop for k from 0 488 for term = 1 then (* term (/ z k)) 489 for sum = (/ ( 1 v)) then (+ sum (/ term (+ k 1 v))) 490 when (< (abs term) (* (abs sum) eps)) 491 return ( (* (gamma ( 1 v)) (expt z ( v 1))) 492 sum)))) 497 (if (and (realp v) 498 (= v (ftruncate v))) 499 ;; v is an integer 500 (let ((n (truncate v))) 501 ( (* (/ (expt z ( v 1)) 502 (gamma v)) 503 ( (psi v) (log z))) 504 (loop for k from 0 below n 505 for term = 1 then (* term (/ z k)) 506 for sum = (/ ( 1 v)) then (+ sum (/ term (+ k 1 v))) 507 when (< (abs term) (* (abs sum) eps)) 508 return sum) 509 (loop for k from n 510 for term = 1 then (* term (/ z k)) 511 for sum = 0 then (+ sum (/ term (+ k 1 v))) 512 when (< (abs term) (* (abs sum) eps)) 513 return sum))) 514 (loop for k from 0 515 for term = 1 then (* term (/ z k)) 516 for sum = (/ ( 1 v)) then (+ sum (/ term (+ k 1 v))) 517 when (< (abs term) (* (abs sum) eps)) 518 return ( (* (gamma ( 1 v)) (expt z ( v 1))) 519 sum))))) 493 520 494 521 (defun expintegrale (v z) … … 680 707 (realpart (ci z)) 681 708 (ci z)))) 709 710 (defconstant bernvalues 711 (makearray 55 712 :initialcontents 713 '(1 714 1/2 715 1/6 716 0 717 1/30 718 0 719 1/42 720 0 721 1/30 722 0 723 5/66 724 0 725 691/2730 726 0 727 7/6 728 0 729 3617/510 730 0 731 43867/798 732 0 733 174611/330 734 0 735 854513/138 736 0 737 236364091/2730 738 0 739 8553103/6 740 0 741 23749461029/870 742 0 743 8615841276005/14322 744 0 745 7709321041217/510 746 0 747 2577687858367/6 748 0 749 26315271553053477373/1919190 750 0 751 2929993913841559/6 752 0 753 261082718496449122051/13530 754 0 755 1520097643918070802691/1806 756 0 757 27833269579301024235023/690 758 0 759 596451111593912163277961/282 760 0 761 5609403368997817686249127547/46410 762 0 763 495057205241079648212477525/66 764 0 765 801165718135489957347924991853/1590 766 0 767 29149963634884862421418123812691/798 768 ))) 769 770 (defun bern (k) 771 (aref bernvalues k)) 772 773 (defun psi (z) 774 "Digamma function defined by 775 776  %gamma + sum(1/k1/(k+z1), k, 1, inf) 777 778 where %gamma is Euler's constant" 779 780 ;; A&S 6.3.7: Reflection formula 781 ;; 782 ;; psi(1z) = psi(z) + %pi*cot(%pi*z) 783 ;; 784 ;; A&S 6.3.6: Recurrence formula 785 ;; 786 ;; psi(n+z) = 1/(z+n1)+1/(z+n2)+...+1/(z+2)+1/(1+z)+psi(1+z) 787 ;; 788 ;; A&S 6.3.8: Asymptotic formula 789 ;; 790 ;; psi(z) ~ log(z)  sum(bern(2*n)/(2*n*z^(2*n)), n, 1, inf) 791 ;; 792 ;; So use reflection formula if Re(z) < 0. For z > 0, use the recurrence 793 ;; formula to increase the argument and then apply the asymptotic formula. 794 795 (cond ((minusp (realpart z)) 796 ( (psi ( 1 z)) 797 (* +pi+ (/ (tan (* +pi+ z)))))) 798 (t 799 (let* ((k (* 2 (1+ (floor (* .41 ( (log (epsilon z) 10))))))) 800 (m 0) 801 (y (expt (+ z k) 2)) 802 (x 0)) 803 (loop for i from 1 upto (floor k 2) do 804 (progn 805 (incf m (+ (/ (+ z i i 1)) 806 (/ (+ z i i 2)))) 807 (setf x (/ (+ x (/ (bern (+ k 2 (* 2 i))) 808 ( k i i 2))) 809 y)))) 810 ( (log (+ z k)) 811 (/ (* 2 (+ z k))) 812 x 813 m))))) 814
Note: See TracChangeset
for help on using the changeset viewer.