Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it

563603 views
Tweet
Share
Share
1
#include "typedef.h"

2
#include "gmp.h"

3
/* #include "gmp-impl.h"*/

4
#include "longtools.h"

5
#include "matrix.h"

6
/**************************************************************************\

7
@---------------------------------------------------------------------------

8
@---------------------------------------------------------------------------

9
@ FILE: long_mat_inv.c

10
@---------------------------------------------------------------------------

11
@---------------------------------------------------------------------------

12
@

13
\**************************************************************************/

14


15
/************************************************************************\

16
@ matrix_TYP *long_mat_inv(A)

17
@ matrix_TYP *A;

18
@

19
@ long_mat_inv(A) calculates the matrix A^{-1} using GNU MP

20
\************************************************************************/

21
matrix_TYP *long_mat_inv(A)

22
matrix_TYP *A;

23
{

24
   MP_INT ***E, **MA, **MB;

25
   MP_INT Ekgv;

26
   int Ecols;

27
   matrix_TYP *erg;

28
   int i,j;

29


30


31
   if(A->cols != A->rows)

32
   {

33
     printf("error: can't invert non-square matrix\n");

34
     exit(3);

35
   }

36
   MA = matrix_to_MP_mat(A);

37
   if((MB = (MP_INT **) malloc(A->cols *sizeof(MP_INT *))) == NULL)

38
   {

39
     printf("malloc of 'MB' in 'long_mat_inv' failed\n");

40
     exit(2);

41
   }

42
   for(i=0;i<A->cols;i++)

43
   {

44
     if((MB[i] = (MP_INT *) malloc(A->cols *sizeof(MP_INT))) == NULL)

45
     {

46
       printf("malloc of 'MB[%d]' in 'long_mat_inv' failed\n", i);

47
       exit(2);

48
     }

49
     for(j=0;j<A->cols;j++)

50
       mpz_init(&MB[i][j]);

51
     mpz_set_si(&MB[i][i], 1);

52
   }

53
   mpz_init(&Ekgv);

54
   E = MP_solve_mat(MA, A->rows, A->cols, MB, A->cols, &Ecols, &Ekgv);

55
   for(i=0;i<A->cols;i++)

56
   {

57
      for(j=0;j<A->cols;j++)

58
      { mpz_clear(&MA[i][j]); mpz_clear(&MB[i][j]);}

59
      free(MA[i]);

60
      free(MB[i]);

61
   }

62
   free(MA);

63
   free(MB);

64


65
   if(E[1] != NULL)

66
   {

67
     printf("Error: matrix in 'long_mat_inv' is singular\n");

68
     exit(3);

69
   }

70
   if(E[0] == NULL)

71
   {

72
     printf("Error in 'long_mat_inv'\n");

73
     exit(3);

74
   }

75
   erg = MP_mat_to_matrix(E[0], A->cols, A->cols);

76
   if(abs(Ekgv._mp_size) > 1)

77
   {

78
     printf("Error: Integer overflow in 'MP_mat_to_matrix'\n");

79
     exit(3);

80
   }

81
   erg->kgv = mpz_get_si(&Ekgv);

82
   for(i=0;i<A->cols;i++)

83
   {

84
      for(j=0;j<A->cols;j++)

85
       mpz_clear(&E[0][i][j]);

86
      free(E[0][i]);

87
   }

88
   free(E[0]);

89
   free(E);

90
   mpz_clear(&Ekgv);

91
   if(A->kgv != 1)

92
   {

93
     for(i=0;i<A->cols;i++)

94
       for(j=0;j<A->cols;j++)

95
         erg->array.SZ[i][j] *= A->kgv;

96
   }

97
   Check_mat(erg);

98
   return(erg);

99
}

100


101