# Ticket #6397: trac_6397-power-series-newton-ncalexan-1.pyx

File trac_6397-power-series-newton-ncalexan-1.pyx, 4.3 KB (added by , 10 years ago) |
---|

Line | |
---|---|

1 | # -*- pyrex -*- |

2 | r""" |

3 | Newton's method for polynomials over power series. |

4 | |

5 | An initial guess `x_0` is refined by the rule `x_{k+1} = x_{k} - |

6 | f(x_k)/f'(x_k)`, the prime denoting the derivative. |

7 | |

8 | AUTHORS: |

9 | |

10 | - Nick Alexander (original implementation) |

11 | """ |

12 | |

13 | # from sage.rings.all import QQ, CDF |

14 | |

15 | from sage.rings.power_series_ring_element cimport PowerSeries |

16 | # from sage.rings.all import infinity |

17 | import sage.misc.misc |

18 | |

19 | def ps_roots(f, prec=5, all=False): |

20 | r""" |

21 | Find roots of polynomials over power series rings like `R[[t]][y]`. |

22 | |

23 | ALGORITHM: |

24 | |

25 | Newton's method. An initial guess `x_0` is refined by the rule `x_{k+1} = |

26 | x_{k} - f(x_k)/f'(x_k)`, the prime denoting the derivative. |

27 | |

28 | WARNING: |

29 | |

30 | .. the algorithm is not known to be numerically stable. Inexact rings |

31 | will suffer from the usual maladies of root finding: see :meth:`roots` for |

32 | more information. |

33 | |

34 | AUTHORS: |

35 | |

36 | - Nick Alexander (original implementation) |

37 | |

38 | EXAMPLES:: |

39 | |

40 | sage: R.<t> = QQ[[]] |

41 | sage: S.<y> = R[] |

42 | |

43 | sage: f = y^3 - (t^2 + t)*y + t^5 + t + 1; g = ps_roots(f, prec=5) |

44 | sage: g |

45 | -1 - 2/3*t - 1/9*t^2 - 1/81*t^3 + 1/27*t^4 + O(t^5) |

46 | sage: f(g) |

47 | O(t^5) |

48 | sage: g = ps_roots(f, prec=12) |

49 | sage: g |

50 | -1 - 2/3*t - 1/9*t^2 - 1/81*t^3 + 1/27*t^4 - 29/81*t^5 + 2284/6561*t^6 - 4624/19683*t^7 + 10693/59049*t^8 - 284333/1594323*t^9 + 50659/177147*t^10 - 2003644/4782969*t^11 + O(t^12) |

51 | sage: f(g) |

52 | O(t^12) |

53 | |

54 | This one should fail. There's an approximation to order 0, namely `y=0 + |

55 | O(t)`, but there is no approximation to order 1. Correspondingly, the |

56 | ``ValuationError`` witnesses that the derivative of `f` at 0 is 0 (to order |

57 | 1) :: |

58 | |

59 | sage: f = y^3 - (t^2 + t)*y + t^5 + t + 11; g = ps_roots(f, prec=5) |

60 | Traceback (most recent call last): |

61 | ... |

62 | ValueError: There exists no approximation to order 0 -- that is, no root of y^3 + 11 defined over Rational Field |

63 | |

64 | With this one, the first order approximation `y=0 + O(t)` cannot be |

65 | improved, since `0 + O(t)` is a root of the derivative of `y^3`, namely `3 |

66 | y^2`:: |

67 | |

68 | sage: f = y^3 - (t^2 + t)*y + t^5 + t; g = ps_roots(f, prec=5) |

69 | Traceback (most recent call last): |

70 | ... |

71 | ValueError: There exists no approximation to order 2 |

72 | |

73 | Now let's test the ``all`` parameter. The first difference is that the |

74 | following doesn't raise an error, it returns an empty list:: |

75 | |

76 | sage: f = y^3 - (2*t^2 + t + 2)^2; ps_roots(f, prec=10, all=True) |

77 | [] |

78 | |

79 | sage: f = y^2 - (2*t^2 + t + 2)^2; gs = ps_roots(f, prec=10, all=True) |

80 | sage: gs |

81 | [2 + t + 2*t^2 + O(t^10), -2 - t - 2*t^2 + O(t^10)] |

82 | sage: [ f(g) for g in gs ] |

83 | [O(t^10), O(t^10)] |

84 | |

85 | sage: f = (y - (2 + t + t^5))*(y - (1 + 15*t + t^2))*(y - (3 + 10*t + t^5)) |

86 | sage: gs = ps_roots(f, prec=10, all=True); gs |

87 | [3 + 10*t + t^5 + O(t^10), 2 + t + t^5 + O(t^10), 1 + 15*t + t^2 + O(t^10)] |

88 | sage: [ f(g) for g in gs ] |

89 | [O(t^10), O(t^10), O(t^10)] |

90 | |

91 | In this example, two roots agree to order 1. That means Newton's method |

92 | doesn't work:: |

93 | |

94 | sage: f = (y - (2 + t + t^5))*(y - (2 + 15*t + t^2))*(y - (3 + 10*t + t^5)) |

95 | sage: ps_roots(f, prec=10, all=True) |

96 | Traceback (most recent call last): |

97 | ... |

98 | ValueError: There exists no approximation to order 2 |

99 | """ |

100 | fp = f.derivative() |

101 | poly_ring, y = f.parent().objgen() |

102 | series_ring, t = poly_ring.base_ring().objgen() |

103 | base_ring = series_ring.base_ring() |

104 | |

105 | # we start with a root of f mod t defined over R |

106 | g = base_ring['y']([ c[0] for c in f ]) # substitute t = 0 |

107 | rts = g.roots(multiplicities=False) |

108 | |

109 | cdef PowerSeries s |

110 | if all is False: |

111 | if not rts: |

112 | raise ValueError, "There exists no approximation to order 0 -- that is, no root of %s defined over %s" % (g, base_ring) |

113 | rts = [ rts[0] ] # just choose a random root |

114 | |

115 | ss = [] |

116 | for rt in rts: |

117 | s = series_ring(rt) |

118 | s._prec = 1 |

119 | |

120 | for cur_prec in sage.misc.misc.newton_method_sizes(prec)[1:]: |

121 | old_prec = s._prec |

122 | s._prec = cur_prec |

123 | D = fp(s) |

124 | if D.valuation() >= old_prec: |

125 | raise ValueError, "There exists no approximation to order %s" % cur_prec |

126 | |

127 | s = s - f(s)/fp(s) |

128 | ss.append(s) |

129 | if all is False: |

130 | return ss[0] |

131 | return ss |