[Pkg-octave-devel] Bug#418158: Bug#418158: more info
Rafael Laboissiere
rafael at debian.org
Sun Feb 24 17:23:30 UTC 2008
* Ólafur Jens Sigurðsson <ojsbug at gmail.com> [2008-02-24 12:34]:
> I just tried some values and the limit seems to be poisscdf(604,604),
> after that the function fails, that is poisscdf(605,605) fails and
> poisscdf(604,604) succeeds.
It fails for smaller values than that, like poisscdf (599, 600). The
culprit is the gammainc function, which ultimately calls the function
defined in d9lgit.f.
I ran the following to discover the "real" limits :
X = [500 : 650];
A = [500 : 650];
nX = length (X);
nA = length (A);
fail = sparse (nX, nA);
for i = 1 : nX
for j = 1 : nA
try
gammainc (X (i), A (j));
catch
fail (i, j) = 1;
end_try_catch
endfor
endfor
[i, j] = spfind (fail);
[(X (i)); (A(j))]'
The result is attached below. It seems that gammainc does not like
arguments X and A close to each other, but this happens only for X and A
starting at 600. For instance, poisscdf (599, 599) and poisscdf (1000,
1010) work fine.
I tried trivial patches as the one below, but this only pushes the limits to
higher values.
--- octave3.0-3.0.0.orig/libcruft/slatec-fn/d9lgit.f
+++ octave3.0-3.0.0/libcruft/slatec-fn/d9lgit.f
@@ -46,7 +46,7 @@
R = 0.D0
P = 1.D0
S = P
- DO 20 K=1,200
+ DO 20 K=1,1000
FK = K
T = (A+FK)*X*(1.D0+R)
R = T/((AX+FK)*(A1X+FK)-T)
@@ -55,7 +55,7 @@
IF (ABS(P).LT.EPS*S) GO TO 30
20 CONTINUE
CALL XERMSG ('SLATEC', 'D9LGIT',
- + 'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 3, 2)
+ + 'NO CONVERGENCE IN 1000 TERMS OF CONTINUED FRACTION', 3, 2)
C
30 HSTAR = 1.0D0 - X*S/A1X
IF (HSTAR .LT. SQEPS) CALL XERMSG ('SLATEC', 'D9LGIT',
--
Rafael
-------------- next part --------------
ans =
600 600
601 601
602 602
603 603
604 604
605 605
605 606
606 606
606 607
607 607
607 608
608 608
608 609
609 609
609 610
610 610
610 611
611 611
610 612
611 612
612 612
611 613
612 613
613 613
612 614
613 614
614 614
613 615
614 615
615 615
614 616
615 616
616 616
615 617
616 617
617 617
615 618
616 618
617 618
618 618
616 619
617 619
618 619
619 619
617 620
618 620
619 620
620 620
618 621
619 621
620 621
621 621
619 622
620 622
621 622
622 622
620 623
621 623
622 623
623 623
620 624
621 624
622 624
623 624
624 624
621 625
622 625
623 625
624 625
625 625
622 626
623 626
624 626
625 626
626 626
623 627
624 627
625 627
626 627
627 627
624 628
625 628
626 628
627 628
628 628
625 629
626 629
627 629
628 629
629 629
626 630
627 630
628 630
629 630
630 630
626 631
627 631
628 631
629 631
630 631
631 631
627 632
628 632
629 632
630 632
631 632
632 632
628 633
629 633
630 633
631 633
632 633
633 633
629 634
630 634
631 634
632 634
633 634
634 634
630 635
631 635
632 635
633 635
634 635
635 635
631 636
632 636
633 636
634 636
635 636
636 636
631 637
632 637
633 637
634 637
635 637
636 637
637 637
632 638
633 638
634 638
635 638
636 638
637 638
638 638
633 639
634 639
635 639
636 639
637 639
638 639
639 639
634 640
635 640
636 640
637 640
638 640
639 640
640 640
635 641
636 641
637 641
638 641
639 641
640 641
641 641
636 642
637 642
638 642
639 642
640 642
641 642
642 642
636 643
637 643
638 643
639 643
640 643
641 643
642 643
643 643
637 644
638 644
639 644
640 644
641 644
642 644
643 644
644 644
638 645
639 645
640 645
641 645
642 645
643 645
644 645
645 645
639 646
640 646
641 646
642 646
643 646
644 646
645 646
646 646
640 647
641 647
642 647
643 647
644 647
645 647
646 647
647 647
641 648
642 648
643 648
644 648
645 648
646 648
647 648
648 648
641 649
642 649
643 649
644 649
645 649
646 649
647 649
648 649
649 649
642 650
643 650
644 650
645 650
646 650
647 650
648 650
649 650
650 650
More information about the Pkg-octave-devel
mailing list