Subject: Extensions and improvements of the Beta funtion
From: Dieter Kaiser
Date: Sun, 08 Feb 2009 00:36:10 +0100
Because I am again working on the Incomplete Beta function I had a look
at the code of the Beta function. These are the points I would like to
suggest to improve:
1. Extension to complex float, bigfloat and complex bigfloat evaluation
Because the Gamma function and the Logarithm of the Gamma function are
fully implemented we can extend the Beta function too.
2. Use log_gamma for numerical evaluation
To avoid an early overflow of an intermediate result use the function
log_gamma and not the function gamma for evaluation.
3. Check the arguments to be valid for numerical evaluation
Check carefully the arguments. The arguments and the sum of the
arguments are not allowed to be a negative integer or a float or
bigfloat representation of a negative integer.
4. Extend expansion for positive integer
There are two changes:
We have one positive integer. We check too if we have a float or
bigfloat representation. With this change we get more consistent
results.
If the second integer is a negative integer, then the sum of the
integers has to be negative to get a valid expansion. This is checked
carefully.
5. Remove code for both arguments a negative integer
The code for the case when both arguments are negative integers don't
work anymore. This code returns a representation of the beta function in
terms of the Factorial function. But the factorial function is not valid
for negative integers.
6. Mirror symmetry and permutation symmetry
I have added mirror symmetry.
I have added permutation symmetry too. But it does not work. I don't no
the reason. Perhaps we have a problem with the noun/verb concept for the
beta function.
Related bug reports are:
SF[831354] beta(-2,1) inconsistent
With the changes we would get as expected:
(%i16) beta(-2,1);
(%o16) -1/2
(%i17) beta(-2.0,1);
(%o17) -0.5
(%i18) beta(-2.0,1.0);
(%o18) -0.5
But this should give -0.5 too (see below):
(%i19) beta(-2,1.0);
(%o19) -1/2
Open problem:
There are some further small inconsistencies for float and bigfloat
numbers e.g. beta(1.0, -2) --> -1/2 and not -0.5. To get all evaluations
very consistent more of the code has to be duplicated for the different
types of numbers.
SF[1860420] "beta_args_sum_to_integer"
As suggested in the bug report it is possible to remove this flag. But I
think the simplification of expressions like beta(1+a,2-a) is
interesting:
(%i21) beta_args_sum_to_integer:true;
(%o21) true
(%i22) beta(1+a,2-a);
(%o22) -%pi*(1-a)*a/(2*sin(%pi*(a+1)))
One possibility would be to allow this simplification generally. But we
have one example in the testsuite which no longer will work.
/* sf bug 1730044 - powerseries(1+x^n,x,0) wrong */
powerseries(1+x^n, x, 0);
('sum(x^(i3*n)/beta(2-i3,i3+1),i3,0,inf))/2;
The result contains a beta function which will simplify to an expression
which generates zeros.
Noun/Verb concept:
The beta function is implemented as a simplifying function, but the
"noun" does not have the name %beta. The noun form is implemented as
$beta. I do not know the reason why the naming convention is changed for
the beta function. But this change causes problems with the parser and
the simplifier.
Perhaps we should change it to the convention %beta for a noun and $beta
for the verb function. This has to be done in about 15 files and at 40
places.
This is the code with the changes:
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Implementation of the Beta function
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defmvar $beta_args_sum_to_integer nil)
;;; The beta funtion has mirror symmetry
(defprop $beta t commutes-with-conjugate)
;;; The Beta function is symmetric beta(a,b) = beta(b,a)
;(eval-when
; #+gcl (load eval)
; #-gcl (:load-toplevel :execute)
; (let (($context '$global) (context '$global))
; (meval '(($declare) $beta $symmetric))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defmfun simpbeta (x vestigial z &aux check)
(declare (ignore vestigial))
(twoargcheck x)
(setq check x)
(let ((u (simpcheck (cadr x) z))
(v (simpcheck (caddr x) z))
($ratprint nil))
(cond ((or (zerop1 u) (zerop1 v))
(if errorsw (throw 'errorsw t) (merror "Zero argument to `beta'")))
;; Check for numerical evaluation in float precision
((and (complex-float-numerical-eval-p u v)
;; We use gamma(u)*gamma(v)/gamma(u+v) for numerical
;; evaluation. Therefore u, v or u+v can not be a
negative
;; integer or a floating point representation of a
negative
;; integer.
(and (or (not (numberp u))
(> u 0)
(not (= (nth-value 1 (truncate u)) 0)))
(and (or (not (numberp v))
(> v 0)
(not (= (nth-value 1 (truncate v)) 0)))
(and (or (not (numberp (add u v)))
(> (add v u) 0)
(not (= (nth-value 1 ($truncate (add u v)))
0)))))))
($rectform (power ($float '$%e)
(add ($log_gamma ($float u))
($log_gamma ($float v))
(mul -1 ($log_gamma ($float (add u v))))))))
;; Check for numerical evaluation in bigfloat precision
((and (not (complex-float-numerical-eval-p u v))
(complex-bigfloat-numerical-eval-p u v)
(and (or (not (mnump u))
(eq ($sign u) '$pos)
(not (eq ($sign (sub ($truncate u) u))
'$zero)))
(or (not (mnump v))
(eq ($sign v) '$pos)
(not (eq ($sign (sub ($truncate v) v))
'$zero)))
(or (not (mnump (add u v)))
(eq ($sign (add u v)) '$pos)
(not (eq ($sign (sub ($truncate (add u v))
(add u v)))
'$zero)))))
($rectform
(power ($bfloat'$%e)
(add ($log_gamma ($bfloat u))
($log_gamma ($bfloat v))
(mul -1 ($log_gamma ($bfloat (add u v))))))))
;; One of the arguments is a positive integer or a float or a
bigfloat
; representation. We try to expand.
((or (and (or (and (integerp u)
(plusp u))
(and (mnump u)
(eq ($sign u) '$pos)
(eq ($sign (sub ($truncate u) u)) '$zero)
(setq u ($truncate u))))
(not (and (mnump v)
(eq ($sign (sub ($truncate v) v)) '$zero)
(eq ($sign v) '$neg)
(eq ($sign (add u v)) '$pos))))
(and (or (and (integerp v)
(plusp v))
(and (mnump v)
(eq ($sign v) '$pos)
(eq ($sign (sub ($truncate v) v)) '$zero)
(setq v ($truncate v))))
(not (and (mnump u)
(eq ($sign (sub ($truncate u) u)) '$zero)
(eq ($sign u) '$neg)
(eq ($sign (add u v)) '$pos)))))
(setq x (add u v))
(power
(mul (sub x 1)
(simplify
(list '(%binomial)
(sub x 2)
(sub (if (and (integerp u) (plusp u)) u v)
1))))
-1))
;; At this point both integers are negative
;; the following code will never work.
; ((and (integerp u) (integerp v))
; (mul2* (div* (list '(mfactorial) (1- u))
; (list '(mfactorial) (+ u v -1)))
; (list '(mfactorial) (1- v))))
;; The sum of the arguments is an integer. We expand.
((or (and (ratnump u)
(ratnump v)
(integerp (setq x (add u v))))
(and $beta_args_sum_to_integer
(integerp (setq x (expand1 (add u v) 1 1)))))
(let ((w (if (symbolp v) v u)))
(div
(mul '$%pi
(simplify (list '(%binomial) (add (1- x) (neg w)) (1-
x))))
(simplify (list '(%sin) (mul w '$%pi))))))
(t (eqtest (list '($beta) u v) check)))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Dieter Kaiser