#
source:
trunk/third/gmp/mpz/gcdext.c
@
18191

Revision 18191, 3.3 KB checked in by ghudson, 22 years ago (diff) |
---|

Line | |
---|---|

1 | /* mpz_gcdext(g, s, t, a, b) -- Set G to gcd(a, b), and S and T such that |

2 | g = as + bt. |

3 | |

4 | Copyright 1991, 1993, 1994, 1995, 1996, 1997, 2000, 2001 Free Software |

5 | Foundation, Inc. |

6 | |

7 | This file is part of the GNU MP Library. |

8 | |

9 | The GNU MP Library is free software; you can redistribute it and/or modify |

10 | it under the terms of the GNU Lesser General Public License as published by |

11 | the Free Software Foundation; either version 2.1 of the License, or (at your |

12 | option) any later version. |

13 | |

14 | The GNU MP Library is distributed in the hope that it will be useful, but |

15 | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |

16 | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |

17 | License for more details. |

18 | |

19 | You should have received a copy of the GNU Lesser General Public License |

20 | along with the GNU MP Library; see the file COPYING.LIB. If not, write to |

21 | the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |

22 | MA 02111-1307, USA. */ |

23 | |

24 | #include <stdio.h> /* for NULL */ |

25 | #include "gmp.h" |

26 | #include "gmp-impl.h" |

27 | |

28 | void |

29 | mpz_gcdext (mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_srcptr a, mpz_srcptr b) |

30 | { |

31 | mp_size_t asize, bsize, usize, vsize; |

32 | mp_srcptr ap, bp; |

33 | mp_ptr up, vp; |

34 | mp_size_t gsize, ssize, tmp_ssize; |

35 | mp_ptr gp, sp, tmp_gp, tmp_sp; |

36 | mpz_srcptr u, v; |

37 | mpz_ptr ss, tt; |

38 | __mpz_struct stmp, gtmp; |

39 | TMP_DECL (marker); |

40 | |

41 | TMP_MARK (marker); |

42 | |

43 | /* mpn_gcdext requires that U >= V. Therefore, we often have to swap U and |

44 | V. This in turn leads to a lot of complications. The computed cofactor |

45 | will be the wrong one, so we have to fix that up at the end. */ |

46 | |

47 | asize = ABS (SIZ (a)); |

48 | bsize = ABS (SIZ (b)); |

49 | ap = PTR (a); |

50 | bp = PTR (b); |

51 | if (asize > bsize || (asize == bsize && mpn_cmp (ap, bp, asize) > 0)) |

52 | { |

53 | usize = asize; |

54 | vsize = bsize; |

55 | up = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); |

56 | vp = (mp_ptr) TMP_ALLOC ((vsize + 1) * BYTES_PER_MP_LIMB); |

57 | MPN_COPY (up, ap, usize); |

58 | MPN_COPY (vp, bp, vsize); |

59 | u = a; |

60 | v = b; |

61 | ss = s; |

62 | tt = t; |

63 | } |

64 | else |

65 | { |

66 | usize = bsize; |

67 | vsize = asize; |

68 | up = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); |

69 | vp = (mp_ptr) TMP_ALLOC ((vsize + 1) * BYTES_PER_MP_LIMB); |

70 | MPN_COPY (up, bp, usize); |

71 | MPN_COPY (vp, ap, vsize); |

72 | u = b; |

73 | v = a; |

74 | ss = t; |

75 | tt = s; |

76 | } |

77 | |

78 | tmp_gp = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); |

79 | tmp_sp = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); |

80 | |

81 | if (vsize == 0) |

82 | { |

83 | tmp_sp[0] = 1; |

84 | tmp_ssize = 1; |

85 | MPN_COPY (tmp_gp, up, usize); |

86 | gsize = usize; |

87 | } |

88 | else |

89 | gsize = mpn_gcdext (tmp_gp, tmp_sp, &tmp_ssize, up, usize, vp, vsize); |

90 | ssize = ABS (tmp_ssize); |

91 | |

92 | PTR (>mp) = tmp_gp; |

93 | SIZ (>mp) = gsize; |

94 | |

95 | PTR (&stmp) = tmp_sp; |

96 | SIZ (&stmp) = (tmp_ssize ^ SIZ (u)) >= 0 ? ssize : -ssize; |

97 | |

98 | if (tt != NULL) |

99 | { |

100 | if (SIZ (v) == 0) |

101 | SIZ (tt) = 0; |

102 | else |

103 | { |

104 | mpz_t x; |

105 | MPZ_TMP_INIT (x, ssize + usize + 1); |

106 | mpz_mul (x, &stmp, u); |

107 | mpz_sub (x, >mp, x); |

108 | mpz_tdiv_q (tt, x, v); |

109 | } |

110 | } |

111 | |

112 | if (ss != NULL) |

113 | { |

114 | if (ALLOC (ss) < ssize) |

115 | _mpz_realloc (ss, ssize); |

116 | sp = PTR (ss); |

117 | MPN_COPY (sp, tmp_sp, ssize); |

118 | SIZ (ss) = SIZ (&stmp); |

119 | } |

120 | |

121 | if (ALLOC (g) < gsize) |

122 | _mpz_realloc (g, gsize); |

123 | gp = PTR (g); |

124 | MPN_COPY (gp, tmp_gp, gsize); |

125 | SIZ (g) = gsize; |

126 | |

127 | TMP_FREE (marker); |

128 | } |

**Note:**See TracBrowser for help on using the repository browser.