sympycore.calculus.functions.elementary

1 # 2 # Created January 2008 by Fredrik Johansson 3 # 4 """ Provides elementary calculus functions sqrt, exp, log, sin, etc and constants pi, E. 5 """ 6 7 __all__ = ['sqrt', 'exp', 'log', 'sin', 'cos', 'tan', 'cot', 'sign', 'E', 'pi', 'gamma'] 8 __docformat__ = "restructuredtext" 9 10 from ..algebra import Calculus, I, NUMBER, TERMS, FACTORS, SYMBOL, TERMS 11 from ..infinity import oo, undefined, CalculusInfinity 12 from ..constants import const_pi, const_E, const_gamma 13 from ..function import Function 14 from ...arithmetic.evalf import evalf 15 from ...arithmetic.numbers import complextypes, realtypes, inttypes 16 from ...arithmetic.number_theory import factorial 17 from ...arithmetic import infinity 18 19 import math 20 21 zero = Calculus.zero 22 one = Calculus.one 23 half = one/2 24 sqrt2 = Calculus.convert('2**(1/2)') 25 sqrt3 = Calculus.convert('3**(1/2)') 26 27 #---------------------------------------------------------------------------# 28 # Exponentials # 29 #---------------------------------------------------------------------------# 30 31 pi = const_pi.as_algebra(Calculus) 32 E = const_E.as_algebra(Calculus) 33 gamma = const_gamma.as_algebra(Calculus) 34 35 Ipi = I*pi 36 Ipi2 = Ipi/2 37

38 -class sign(Function):

39 - def __new__(cls, arg):

40 if not isinstance(arg, Calculus): 41 arg = Calculus.convert(arg) 42 return Calculus(cls, arg)

43

44 -class sqrt(Function):

45 - def __new__(cls, arg):

46 return arg ** half

47

48 -class exp(Function):

49 - def __new__(cls, arg):

50 return E ** arg

51 52 @staticmethod

53 - def derivative(arg):

54 return E ** arg

55 56 log_number_table = { 57 zero.data : -oo, 58 one.data : zero, 59 I.data : Ipi2, 60 (-I).data : -Ipi2 61 } 62

63 -class log(Function):

64 - def __new__(cls, arg, base=E):

65 if type(arg) is not Calculus: 66 if isinstance(arg, CalculusInfinity): 67 if arg == oo: 68 return oo 69 if arg == undefined: 70 return undefined 71 return Calculus(cls, arg) 72 else: 73 arg = Calculus.convert(arg) 74 if base != E: 75 base = Calculus.convert(base) 76 bd = base.data 77 ad = arg.data 78 if base.head is NUMBER and isinstance(bd, inttypes) and \ 79 arg.head is NUMBER and isinstance(ad, inttypes) and \ 80 ad > 0 and bd > 1: 81 l = int(math.log(ad, bd) + 0.5) 82 if bd**l == ad: 83 return Calculus.convert(l) 84 return cls(arg) / cls(base) 85 head = arg.head 86 data = arg.data 87 if head is NUMBER: 88 v = log_number_table.get(data) 89 if v is not None: 90 return v 91 if isinstance(data, realtypes) and data < 0: 92 return Ipi + log(-arg) 93 if isinstance(data, complextypes) and data.real == 0: 94 im = data.imag 95 if im > 0: return Calculus(cls, Calculus(NUMBER, im)) + Ipi2 96 if im < 0: return Calculus(cls, Calculus(NUMBER, -im)) - Ipi2 97 return Calculus(cls, arg) 98 if arg == E: 99 return one 100 from ..relational import is_positive 101 if head is FACTORS and len(data) == 1: 102 base, expt = data.items()[0] 103 if is_positive(base) and isinstance(expt, realtypes): 104 return Calculus(cls, base) * expt 105 if head is TERMS and len(data) == 1: 106 term, coeff = data.items()[0] 107 if (isinstance(coeff, realtypes) and coeff < 0) and is_positive(base): 108 return Ipi + log(-arg) 109 return Calculus(cls, arg)

110 111 @classmethod

112 - def derivative(cls, arg):

113 return 1/arg

114 115 @classmethod

116 - def nth_derivative(cls, arg, n=1):

117 return (-1)**(n-1) * factorial(n-1) * arg**(-n)

118 119 #---------------------------------------------------------------------------# 120 # Trigonometric functions # 121 #---------------------------------------------------------------------------# 122 123 # Tabulated values of sin(x) at multiples of pi/12 124 C0 = (sqrt3-1)/(2*sqrt2) 125 C1 = one/2 126 C2 = sqrt2/2 127 C3 = sqrt3/2 128 C4 = (sqrt3+1)/(2*sqrt2) 129 130 # Replace entries with None to prevent from evaluating 131 sine_table = [ \ 132 Calculus.zero, C0, C1, C2, C3, C4, Calculus.one, C4, C3, C2, C1,C0, 133 Calculus.zero,-C0,-C1,-C2,-C3,-C4,-Calculus.one,-C4,-C3,-C2,-C1,-C0] 134

135 -def get_pi_shift(arg, N):

136 """Parse as x, n where arg = x + n*pi/N""" 137 if arg.head is TERMS: 138 # Optimizing for len == 1 gives 2x speedup in simple cases 139 if len(arg.data) == 1: 140 e, c = arg.data.items()[0] 141 if e == pi: 142 c *= N 143 if isinstance(c, (int, long)): 144 return zero, c 145 c = arg.data.get(pi) 146 if c is not None: 147 c *= N 148 if isinstance(c, (int, long)): 149 d = arg.data.copy() 150 d.pop(pi) 151 return Calculus.Terms(*d.items()), c 152 return arg, 0 153 if arg == pi: 154 return zero, N 155 return arg, 0

156

157 -def has_leading_sign(arg):

158 """Detect symmetry cases for odd/even functions.""" 159 if arg.head is NUMBER: 160 if arg < 0: 161 return True 162 if arg.head is TERMS and len(arg.data) == 1: 163 e, c = arg.data.items()[0] 164 if c < 0: 165 return True 166 return None

167

168 -class TrigonometricFunction(Function):

169 170 parity = None # 'even' or 'odd' 171 period = None # multiple of pi 172

173 - def __new__(cls, arg):

174 if type(arg) is not Calculus: 175 if isinstance(arg, CalculusInfinity): 176 if arg == undefined: 177 return undefined 178 return Calculus(cls, arg) 179 else: 180 arg = Calculus.convert(arg) 181 x, m = get_pi_shift(arg, 12) 182 m %= (12*cls.period) 183 if x == zero: 184 v = cls.eval_direct(m) 185 if v is not None: 186 return v 187 period = cls.period 188 negate_result = False 189 # Full-period symmetry 190 if not m % (12*period): 191 arg = x 192 else: 193 # Half-period symmetry 194 if not m % 12: 195 arg = x 196 negate_result = True 197 # Quarter-period symmetry 198 elif not m % 6: 199 f = conjugates[cls] 200 sign = (-1)**((((m-6)//12) % period) + (f.parity == 'odd')) 201 return sign * f(x) 202 if has_leading_sign(arg): 203 arg = -arg 204 negate_result ^= (cls.parity == 'odd') 205 if negate_result: 206 return -Calculus(cls, arg) 207 else: 208 return Calculus(cls, arg)

209

210 -class sin(TrigonometricFunction):

211 parity = 'odd' 212 period = 2 213 @classmethod

214 - def eval_direct(cls, m):

215 return sine_table[m % 24]

216 217 @classmethod

218 - def derivative(cls, arg, n=1):

219 if n == 1: 220 return cos(arg)

221 222 @classmethod

223 - def nth_derivative(cls, arg, n):

224 return sin(arg + n*pi/2)

225

226 -class cos(TrigonometricFunction):

227 parity = 'even' 228 period = 2 229 230 @classmethod

231 - def eval_direct(cls, m):

232 return sine_table[(m+6) % 24]

233 234 @classmethod

235 - def derivative(cls, arg, n=1):

236 return -sin(arg)

237 238 @classmethod

239 - def nth_derivative(cls, arg, n=1):

240 return cos(arg + n*pi/2)

241

242 -class tan(TrigonometricFunction):

243 parity = 'odd' 244 period = 1 245 246 @classmethod

247 - def eval_direct(cls, m):

248 a = sine_table[m % 24] 249 b = sine_table[(m+6) % 24] 250 if a == None or b == None: 251 return 252 return a / b

253 254 @classmethod

255 - def derivative(cls, arg):

256 return 1+tan(arg)**2

257

258 -class cot(TrigonometricFunction):

259 parity = 'odd' 260 period = 1 261 262 @classmethod

263 - def eval_direct(cls, m):

264 a = sine_table[m % 24] 265 b = sine_table[(m+6) % 24] 266 if a == None or b == None: 267 return 268 return b / a

269 270 @classmethod

271 - def derivative(cls, arg, n=1):

272 return -(1+cot(arg)**2)

273 274 # pi/2-x symmetry 275 conjugates = { 276 sin : cos, 277 cos : sin, 278 tan : cot, 279 cot : tan 280 } 281

Source Code for Module sympycore.calculus.functions.elementary