267 {
329
330
332 if( n==0 )
333 {
335 }
336 if( n==1 )
337 {
338 if( d(1)<0 )
339 {
340 d(1) = -d(1);
342 {
344 }
345 }
347 }
348
349
350
351
352 work0.setbounds(1, n-1);
353 work1.setbounds(1, n-1);
354 work2.setbounds(1, n-1);
355 work3.setbounds(1, n-1);
365
366
367
368
369 etemp.setbounds(1, n);
370 for(
i=1;
i<=n-1;
i++)
371 {
373 }
375 for(
i=1;
i<=n-1;
i++)
376 {
378 }
379 e(n) = 0;
381
382
383
384
387
388
389
390
391
393 {
394 for(
i=1;
i<=n-1;
i++)
395 {
396 rotations::generaterotation<Precision>(d(
i), e(
i),
cs,
sn, r);
402 }
403
404
405
406
408 {
410 }
412 {
414 }
415 }
416
417
418
419
420
421
425 {
427 }
428
429
430
431
434 {
435 smax = amp::maximum<Precision>(
smax, amp::abs<Precision>(d(
i)));
436 }
437 for(
i=1;
i<=n-1;
i++)
438 {
439 smax = amp::maximum<Precision>(
smax, amp::abs<Precision>(e(
i)));
440 }
443 {
444
445
446
447
448 sminoa = amp::abs<Precision>(d(1));
450 {
453 {
454 mu = amp::abs<Precision>(d(
i))*(
mu/(
mu+amp::abs<Precision>(e(
i-1))));
457 {
458 break;
459 }
460 }
461 }
464 }
465 else
466 {
467
468
469
470
472 }
473
474
475
476
477
478
483
484
485
486
488
489
490
491
492 while( true )
493 {
494
495
496
497
499 {
500 break;
501 }
503 {
506 }
507
508
509
510
511 if(
tol<0 && amp::abs<Precision>(d(
m))<=
thresh )
512 {
514 }
515 smax = amp::abs<Precision>(d(
m));
519 {
521 abss = amp::abs<Precision>(d(
ll));
522 abse = amp::abs<Precision>(e(
ll));
524 {
526 }
528 {
530 break;
531 }
534 }
536 {
538 }
539 else
540 {
541
542
543
544
547 {
548
549
550
551
553 continue;
554 }
555 }
557
558
559
560
562 {
563
564
565
566
571
572
573
574
576 {
584 }
586 {
594 }
596 {
604 }
606 continue;
607 }
608
609
610
611
612
613
614
615
616
617
620 {
622 }
624 {
626 }
628 {
629 if( amp::abs<Precision>(d(
ll))>=amp::abs<Precision>(d(
m)) )
630 {
631
632
633
634
636 }
637 else
638 {
639
640
641
642
644 }
645 }
646
647
648
649
651 {
652
653
654
655
656
657 if( amp::abs<Precision>(e(
m-1))<=amp::abs<Precision>(
tol)*amp::abs<Precision>(d(
m)) ||
tol<0 && amp::abs<Precision>(e(
m-1))<=
thresh )
658 {
660 continue;
661 }
663 {
664
665
666
667
668
669 mu = amp::abs<Precision>(d(
ll));
673 {
674 if( amp::abs<Precision>(e(
lll))<=
tol*
mu )
675 {
678 break;
679 }
681 mu = amp::abs<Precision>(d(
lll+1))*(
mu/(
mu+amp::abs<Precision>(e(
lll))));
683 }
685 {
686 continue;
687 }
688 }
689 }
690 else
691 {
692
693
694
695
696
697 if( amp::abs<Precision>(e(
ll))<=amp::abs<Precision>(
tol)*amp::abs<Precision>(d(
ll)) ||
tol<0 && amp::abs<Precision>(e(
ll))<=
thresh )
698 {
700 continue;
701 }
703 {
704
705
706
707
708
709 mu = amp::abs<Precision>(d(
m));
713 {
714 if( amp::abs<Precision>(e(
lll))<=
tol*
mu )
715 {
718 break;
719 }
721 mu = amp::abs<Precision>(d(
lll))*(
mu/(
mu+amp::abs<Precision>(e(
lll))));
723 }
725 {
726 continue;
727 }
728 }
729 }
732
733
734
735
736
738 {
739
740
741
742
743 shift = 0;
744 }
745 else
746 {
747
748
749
750
752 {
753 sll = amp::abs<Precision>(d(
ll));
755 }
756 else
757 {
758 sll = amp::abs<Precision>(d(
m));
760 }
761
762
763
764
766 {
767 if( amp::sqr<Precision>(shift/
sll)<
eps )
768 {
769 shift = 0;
770 }
771 }
772 }
773
774
775
776
778
779
780
781
782 if( shift==0 )
783 {
785 {
786
787
788
789
790
794 {
795 rotations::generaterotation<Precision>(d(
i)*
cs, e(
i),
cs,
sn, r);
797 {
799 }
806 }
810
811
812
813
815 {
817 }
819 {
821 }
823 {
825 }
826
827
828
829
830 if( amp::abs<Precision>(e(
m-1))<=
thresh )
831 {
833 }
834 }
835 else
836 {
837
838
839
840
841
845 {
846 rotations::generaterotation<Precision>(d(
i)*
cs, e(
i-1),
cs,
sn, r);
848 {
850 }
857 }
861
862
863
864
866 {
868 }
870 {
872 }
874 {
876 }
877
878
879
880
881 if( amp::abs<Precision>(e(
ll))<=
thresh )
882 {
884 }
885 }
886 }
887 else
888 {
889
890
891
892
894 {
895
896
897
898
899
903 {
904 rotations::generaterotation<Precision>(
f,
g,
cosr,
sinr, r);
906 {
908 }
913 rotations::generaterotation<Precision>(
f,
g,
cosl,
sinl, r);
918 {
921 }
926 }
928
929
930
931
933 {
935 }
937 {
939 }
941 {
943 }
944
945
946
947
948 if( amp::abs<Precision>(e(
m-1))<=
thresh )
949 {
951 }
952 }
953 else
954 {
955
956
957
958
959
963 {
964 rotations::generaterotation<Precision>(
f,
g,
cosr,
sinr, r);
966 {
968 }
973 rotations::generaterotation<Precision>(
f,
g,
cosl,
sinl, r);
978 {
981 }
986 }
988
989
990
991
992 if( amp::abs<Precision>(e(
ll))<=
thresh )
993 {
995 }
996
997
998
999
1001 {
1003 }
1005 {
1007 }
1009 {
1011 }
1012 }
1013 }
1014
1015
1016
1017
1018 continue;
1019 }
1020
1021
1022
1023
1025 {
1027 {
1029
1030
1031
1032
1034 {
1036 }
1037 }
1038 }
1039
1040
1041
1042
1043
1044 for(
i=1;
i<=n-1;
i++)
1045 {
1046
1047
1048
1049
1052 for(
j=2;
j<=n+1-
i;
j++)
1053 {
1055 {
1058 }
1059 }
1061 {
1062
1063
1064
1065
1069 {
1074 }
1076 {
1081 }
1083 {
1088 }
1089 }
1090 }
1092 }
static const ampf getAlgoPascalEpsilon()
static const ampf getAlgoPascalMinNumber()
void setbounds(int iLow, int iHigh)
raw_vector< T > getrow(int iRow, int iColumnStart, int iColumnEnd)
raw_vector< T > getcolumn(int iColumn, int iRowStart, int iRowEnd)
static matrix mu(matrix A, const ring R)
int maxint(int m1, int m2)
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
void vadd(raw_vector< T > vdst, const_raw_vector< T > vsrc)
void vmul(raw_vector< T > vdst, T2 alpha)
void vsub(raw_vector< T > vdst, const_raw_vector< T > vsrc)