1 r"""

2 Fast Numerical Evaluation.

3

4 For many applications such as numerical integration, differential

5 equation approximation, plotting a 3d surface, optimization problems,

6 monte-carlo simulations, etc., one wishes to pass around and evaluate

7 a single algebraic expression many, many times at various floating

8 point values. Doing this via recursive calls over a python

9 representation of the object (even if maxima or other outside packages

10 are not involved) is extremely inefficient.

11

12 Up until now the solution has been to use lambda expressions, but this

13 is neither intuitive, Sage-like, nor efficient (compared to operating

14 on raw C doubles). This module provides a representation of algebraic

15 expression in Reverse Polish Notation, and provides an efficient

16 interpreter on C double values as a callable python object. It does

17 what it can in C, and will call out to Python if necessary.

18

19 Essential to the understanding of this class is the distinction

20 between symbolic expressions and callable symbolic expressions (where

21 the later binds argument names to argument positions). The

22 \code{*vars} parameter passed around encapsulates this information.

23

24 See the function \code{fast_float(f, *vars)} to create a fast-callable

25 version of f.

26

27 To provide this interface for a class, implement

28 \code{_fast_float_(self, *vars)}. The basic building blocks are

29 provided by the functions \code{fast_float_constant} (returns a

30 constant function), \code{fast_float_arg} (selects the $n$-th value

31 when called with $\ge n$ arguments), and \code{fast_float_func} which

32 wraps a callable Python function. These may be combined with the

33 standard Python arithmetic operators, and support many of the basic

34 math functions such sqrt, exp, and trig functions.

35

36 EXAMPLES:

37 sage: from sage.ext.fast_eval import fast_float

38 sage: f = fast_float(sqrt(x^2+1), 'x')

39 sage: f(1)

40 1.4142135623730951

41 sage: f.op_list()

42 ['load 0', 'push 2.0', 'call pow(2)', 'push 1.0', 'add', 'call sqrt(1)']

43

44 To interpret that last line, we load argument 0 ('x' in this case) onto

45 the stack, push the constant 2.0 onto the stack, call the pow function

46 (which takes 2 arguments from the stack), push the constant 1.0, add the

47 top two arguments of the stack, and then call sqrt.

48

49 Here we take sin of the first argument and add it to f:

50 sage: from sage.ext.fast_eval import fast_float_arg

51 sage: g = fast_float_arg(0).sin()

52 sage: (f+g).op_list()

53 ['load 0', 'push 2.0', 'call pow(2)', 'push 1.0', 'add', 'call sqrt(1)', 'load 0', 'call sin(1)', 'add']

54

55 AUTHOR:

56 -- Robert Bradshaw (2008-10): Initial version

57 """

58

59

60 #*****************************************************************************

61 # Copyright (C) 2008 Robert Bradshaw <robertwb@math.washington.edu>

62 #

63 # Distributed under the terms of the GNU General Public License (GPL)

64 #

65 # This code is distributed in the hope that it will be useful,

66 # but WITHOUT ANY WARRANTY; without even the implied warranty of

67 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

68 # General Public License for more details.

69 #

70 # The full text of the GPL is available at:

71 #

72 # http://www.gnu.org/licenses/

73 #*****************************************************************************

74

75 include "stdsage.pxi"

76

77 cdef extern from "Python.h":

78 int PyInt_AS_LONG(PyObject*)

79 PyObject* PyTuple_New(Py_ssize_t size)

80 PyObject* PyTuple_GET_ITEM(PyObject* t, Py_ssize_t index)

81 void PyTuple_SET_ITEM(PyObject* t, Py_ssize_t index, PyObject* item)

82 object PyObject_CallObject(PyObject* func, PyObject* args)

83 PyObject* PyFloat_FromDouble(double d)

84 void Py_DECREF(PyObject *)

85

86 cdef extern from "math.h":

87 double sqrt(double)

88 double pow(double, double)

89

90 double ceil(double)

91 double floor(double)

92

93 double sin(double)

94 double cos(double)

95 double tan(double)

96

97 double asin(double)

98 double acos(double)

99 double atan(double)

100 double atan2(double, double)

101

102 double sinh(double)

103 double cosh(double)

104 double tanh(double)

105

106 double asinh(double)

107 double acosh(double)

108 double atanh(double)

109

110 double exp(double)

111 double log(double)

112 double log2(double)

113 double log10(double)

114

115 cdef extern from *:

116 void* memcpy(void* dst, void* src, size_t len)

117

118 cdef inline int max(int a, int b):

119 return a if a > b else b

120

121 cdef inline int min(int a, int b):

122 return a if a < b else b

123

124 cdef enum:

125 # stack

126 LOAD_ARG # push input argument n onto the stack

127 PUSH_CONST

128 POP

129 POP_N

130 DUP

131

132 # basic arithamtic

133 ADD

134 SUB

135 MUL

136 DIV

137 NEG

138 ABS

139 INVERT

140

141 # functional

142 ONE_ARG_FUNC

143 TWO_ARG_FUNC

144 PY_FUNC

145

146

147 # These two dictionaries are for printable and machine independant representation.

148

149 op_names = {

150 LOAD_ARG: 'load',

151 PUSH_CONST: 'push',

152 POP: 'pop',

153 POP_N: 'popn',

154 DUP: 'dup',

155

156 ADD: 'add',

157 SUB: 'sub',

158 MUL: 'mul',

159 DIV: 'div',

160 NEG: 'neg',

161 ABS: 'abs',

162 INVERT: 'invert',

163

164 ONE_ARG_FUNC: 'call',

165 TWO_ARG_FUNC: 'call',

166 PY_FUNC: 'py_call',

167 }

168

169 cfunc_names = {

170 <size_t>&sqrt: 'sqrt',

171 <size_t>&pow: 'pow',

172

173 <size_t>&ceil: 'ceil',

174 <size_t>&floor: 'floor',

175

176 <size_t>&sin: 'sin',

177 <size_t>&cos: 'cos',

178 <size_t>&tan: 'tan',

179

180 <size_t>&asin: 'asin',

181 <size_t>&atan: 'atan',

182 <size_t>&atan2: 'atan2',

183

184 <size_t>&sinh: 'sinh',

185 <size_t>&cosh: 'cosh',

186 <size_t>&tanh: 'tanh',

187

188 <size_t>&asinh: 'asinh',

189 <size_t>&acosh: 'acosh',

190 <size_t>&atanh: 'atanh',

191

192 <size_t>&exp: 'exp',

193 <size_t>&log: 'log',

194 <size_t>&log2: 'log2',

195 <size_t>&log10: 'log10',

196

197 }

198

199 cdef reverse_map(m):

200 r = {}

201 for key, value in m.iteritems():

202 r[value] = key

203 return r

204

205 # With all the functionality around the op struct, perhaps there should be

206 # a wrapper class, though we still wish to operate on pure structs for speed.

207

208 cdef op_to_string(fast_double_op op):

209 s = op_names[op.type]

210 if op.type in [LOAD_ARG, POP_N]:

211 s += " %s" % op.params.n

212 elif op.type == PUSH_CONST:

213 s += " %s" % op.params.c

214 elif op.type in [ONE_ARG_FUNC, TWO_ARG_FUNC]:

215 try:

216 cname = cfunc_names[<size_t>op.params.func]

217 except KeyError:

218 cname = "0x%x" % <size_t>op.params.func

219 s += " %s(%s)" % (cname, 1 if op.type == ONE_ARG_FUNC else 2)

220 elif op.type == PY_FUNC:

221 n, func = <object>(op.params.func)

222 s += " %s(%s)" % (func, n)

223 return s

224

225 cdef op_to_tuple(fast_double_op op):

226 s = op_names[op.type]

227 if op.type in [LOAD_ARG, POP_N]:

228 param = op.params.n

229 elif op.type == PUSH_CONST:

230 param = op.params.c

231 elif op.type in [ONE_ARG_FUNC, TWO_ARG_FUNC]:

232 param_count = 1 if op.type == ONE_ARG_FUNC else 2

233 try:

234 param = param_count, cfunc_names[<size_t>op.params.func]

235 except KeyError:

236 raise ValueError, "Unknown C function: 0x%x" % <size_t>op.params.func

237 elif op.type == PY_FUNC:

238 param = <object>(op.params.func)

239 else:

240 param = None

241 if param is None:

242 return (s,)

243 else:

244 return s, param

245

246 def _unpickle_FastDoubleFunc(nargs, max_height, op_list):

247 cdef FastDoubleFunc self = PY_NEW(FastDoubleFunc)

248 self.nops = len(op_list)

249 self.nargs = nargs

250 self.max_height = max_height

251 self.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op) * self.nops)

252 self.allocate_stack()

253 cfunc_addresses = reverse_map(cfunc_names)

254 op_enums = reverse_map(op_names)

255 cdef size_t address

256 cdef int i = 0, type

257 for op in op_list:

258 self.ops[i].type = type = op_enums[op[0]]

259 if type in [LOAD_ARG, POP_N]:

260 self.ops[i].params.n = op[1]

261 elif type == PUSH_CONST:

262 self.ops[i].params.c = op[1]

263 elif type in [ONE_ARG_FUNC, TWO_ARG_FUNC]:

264 param_count, cfunc = op[1]

265 address = cfunc_addresses[cfunc]

266 self.ops[i].params.func = <void *>address

267 self.ops[i].type = ['', ONE_ARG_FUNC, TWO_ARG_FUNC][param_count]

268 elif type == PY_FUNC:

269 if self.py_funcs is None:

270 self.py_funcs = op[1]

271 else:

272 self.py_funcs = self.py_funcs + (op[1],)

273 self.ops[i].params.func = <void *>op[1]

274 i += 1

275 return self

276

277

278 # This is where we wish we had case statements...

279 # It looks like gcc might be smart enough to figure it out.

280 cdef inline int process_op(fast_double_op op, double* stack, double* argv, int top) except -2:

281

282 # print [stack[i] for i from 0 <= i <= top], ':', op.type

283

284 cdef int i, n

285 cdef PyObject* py_args

286

287 # We have to do some trickery because Pyrex dissallows function pointer casts

288 # This will be removed in a future version of Cython.

289 cdef double (*f)(double)

290 cdef void** fp = <void **>&f

291 cdef double (*ff)(double, double)

292 cdef void** ffp = <void **>&ff

293

294 if op.type == LOAD_ARG:

295 stack[top+1] = argv[op.params.n]

296 return top+1

297

298 elif op.type == PUSH_CONST:

299 stack[top+1] = op.params.c

300 return top+1

301

302 elif op.type == POP:

303 return top-1

304

305 elif op.type == POP_N:

306 return top-op.params.n

307

308 elif op.type == DUP:

309 stack[top+1] = stack[top]

310 return top+1

311

312 elif op.type == ADD:

313 stack[top-1] += stack[top]

314 return top-1

315

316 elif op.type == SUB:

317 stack[top-1] -= stack[top]

318 return top-1

319

320 elif op.type == MUL:

321 stack[top-1] *= stack[top]

322 return top-1

323

324 elif op.type == DIV:

325 stack[top-1] /= stack[top]

326 return top-1

327

328 elif op.type == NEG:

329 stack[top] = -stack[top]

330 return top

331

332 elif op.type == ABS:

333 if stack[top] < 0:

334 stack[top] = -stack[top]

335 return top

336

337 elif op.type == INVERT:

338 stack[top] = 1/stack[top]

339 return top

340

341 elif op.type == ONE_ARG_FUNC:

342 fp[0] = op.params.func

343 stack[top] = f(stack[top])

344 return top

345

346 elif op.type == TWO_ARG_FUNC:

347 ffp[0] = op.params.func

348 stack[top-1] = ff(stack[top-1], stack[top])

349 return top-1

350

351 elif op.type == PY_FUNC:

352 # Even though it's python, optimize this because it'll be used often...

353 # We also want to avoid cluttering the header and footer

354 # of this function Python variables bring.

355 n = PyInt_AS_LONG(PyTuple_GET_ITEM(op.params.func, 0))

356 top = top - n + 1

357 py_args = PyTuple_New(n)

358 for i from 0 <= i < n:

359 PyTuple_SET_ITEM(py_args, i, PyFloat_FromDouble(stack[top+i]))

360 stack[top] = PyObject_CallObject(PyTuple_GET_ITEM(op.params.func, 1), py_args)

361 Py_DECREF(py_args)

362 return top

363

364 raise RuntimeError, "Bad op code %s" % op.type

365

366

367 cdef class FastDoubleFunc:

368 """

369 This class is for fast evaluation of algebraic expressions over

370 the real numbers (e.g. for plotting). It represents an expression

371 as a stack-based series of operations.

372

373 EXAMPLES:

374 sage: from sage.ext.fast_eval import FastDoubleFunc

375 sage: f = FastDoubleFunc('const', 1.5) # the constant function

376 sage: f()

377 1.5

378 sage: g = FastDoubleFunc('arg', 0) # the first argument

379 sage: g(5)

380 5.0

381 sage: h = f+g

382 sage: h(17)

383 18.5

384 sage: h = h.sin()

385 sage: h(pi/2-1.5)

386 1.0

387 sage: h.is_pure_c()

388 True

389 sage: list(h)

390 ['push 1.5', 'load 0', 'add', 'call sin(1)']

391

392 We can wrap Python functions too:

393 sage: h = FastDoubleFunc('callable', lambda x,y: x*x*x - y, g, f)

394 sage: h(10)

395 998.5

396 sage: h.is_pure_c()

397 False

398 sage: list(h) # random address

399 ['load 0', 'push 1.5', 'py_call <function <lambda> at 0x9fedf70>(2)']

400

401 Here's a more complicated expression:

402 sage: from sage.ext.fast_eval import fast_float_constant, fast_float_arg

403 sage: a = fast_float_constant(1.5)

404 sage: b = fast_float_constant(3.14)

405 sage: c = fast_float_constant(7)

406 sage: x = fast_float_arg(0)

407 sage: y = fast_float_arg(1)

408 sage: f = a*x^2 + b*x + c - y/sqrt(sin(y)^2+a)

409 sage: f(2,3)

410 16.846610528508116

411 sage: f.max_height

412 4

413 sage: f.is_pure_c()

414 True

415 sage: list(f)

416 ['push 1.5', 'load 0', 'dup', 'mul', 'mul', 'push 3.14', 'load 0', 'mul', 'add', 'push 7.0', 'add', 'load 1', 'load 1', 'call sin(1)', 'dup', 'mul', 'push 1.5', 'add', 'call sqrt(1)', 'div', 'sub']

417

418

419 AUTHOR:

420 -- Robert Bradshaw

421 """

422 def __init__(self, type, param, *args):

423

424 cdef FastDoubleFunc arg

425 cdef int i

426

427 if type == 'arg':

428 self.nargs = param+1

429 self.nops = 1

430 self.max_height = 1

431 self.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op))

432 self.ops[0].type = LOAD_ARG

433 self.ops[0].params.n = param

434

435 elif type == 'const':

436 self.nargs = 0

437 self.nops = 1

438 self.max_height = 1

439 self.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op))

440 self.ops[0].type = PUSH_CONST

441 self.ops[0].params.c = param

442

443 elif type == 'callable':

444 py_func = len(args), param

445 self.py_funcs = (py_func,) # just so it doesn't get garbage collected

446 self.nops = 1

447 self.nargs = 0

448 for i from 0 <= i < len(args):

449 a = args[i]

450 if not isinstance(a, FastDoubleFunc):

451 a = FastDoubleFunc('const', a)

452 args = args[:i] + (a,) + args[i+1:]

453 arg = a

454 self.nops += arg.nops

455 if arg.py_funcs is not None:

456 self.py_funcs += arg.py_funcs

457 self.nargs = max(self.nargs, arg.nargs)

458 self.max_height = max(self.max_height, arg.max_height+i)

459 self.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op) * self.nops)

460 if self.ops == NULL:

461 raise MemoryError

462 i = 0

463 for arg in args:

464 memcpy(self.ops + i, arg.ops, sizeof(fast_double_op) * arg.nops)

465 i += arg.nops

466 self.ops[self.nops-1].type = PY_FUNC

467 self.ops[self.nops-1].params.func = <void *>py_func

468

469 else:

470 raise ValueError, "Unknown operation: %s" % type

471

472 self.allocate_stack()

473

474

475 cdef int allocate_stack(FastDoubleFunc self) except -1:

476 self.argv = <double*>sage_malloc(sizeof(double) * self.nargs)

477 if self.argv == NULL:

478 raise MemoryError

479 self.stack = <double*>sage_malloc(sizeof(double) * self.max_height)

480 if self.stack == NULL:

481 raise MemoryError

482

483 def __del__(self):

484 if self.ops:

485 sage_free(self.ops)

486 if self.stack:

487 sage_free(self.stack)

488 if self.argv:

489 sage_free(self.argv)

490

491 def __reduce__(self):

492 """

493 TESTS:

494 sage: from sage.ext.fast_eval import fast_float_arg, fast_float_func

495 sage: f = fast_float_arg(0).sin() * 10 + fast_float_func(hash, fast_float_arg(1))

496 sage: loads(dumps(f)) == f

497 True

498 """

499 L = [op_to_tuple(self.ops[i]) for i from 0 <= i < self.nops]

500 return _unpickle_FastDoubleFunc, (self.nargs, self.max_height, L)

501

502 def __cmp__(self, other):

503 """

504 Two functions are considered equal if they represent the same

505 exact sequence of operations.

506

507 TESTS:

508 sage: from sage.ext.fast_eval import fast_float_arg

509 sage: fast_float_arg(0) == fast_float_arg(0)

510 True

511 sage: fast_float_arg(0) == fast_float_arg(1)

512 False

513 sage: fast_float_arg(0) == fast_float_arg(0).sin()

514 False

515 """

516 cdef int c, i

517 cdef FastDoubleFunc left, right

518 try:

519 left, right = self, other

520 c = cmp((left.nargs, left.nops, left.max_height),

521 (right.nargs, right.nops, right.max_height))

522 if c != 0:

523 return c

524 for i from 0 <= i < self.nops:

525 if left.ops[i].type != right.ops[i].type:

526 return cmp(left.ops[i].type, right.ops[i].type)

527 for i from 0 <= i < self.nops:

528 c = cmp(op_to_tuple(left.ops[i]), op_to_tuple(right.ops[i]))

529 if c != 0:

530 return c

531 return c

532 except TypeError:

533 return cmp(type(self), type(other))

534

535 def __call__(FastDoubleFunc self, *args):

536 """

537 EXAMPLES:

538 sage: from sage.ext.fast_eval import fast_float_arg

539 sage: f = fast_float_arg(2)

540 sage: f(0,1,2,3)

541 2.0

542 sage: f(10)

543 Traceback (most recent call last):

544 ...

545 TypeError: Wrong number of arguments (need at least 3, got 1)

546 sage: f('blah', 1, 2, 3)

547 Traceback (most recent call last):

548 ...

549 TypeError: a float is required

550 """

551 if len(args) < self.nargs:

552 raise TypeError, "Wrong number of arguments (need at least %s, got %s)" % (self.nargs, len(args))

553 cdef int i = 0

554 for i from 0 <= i < self.nargs:

555 self.argv[i] = args[i]

556 res = self._call_c(self.argv)

557 return res

558

559 cdef double _call_c(FastDoubleFunc self, double* argv) except? -2:

560 # The caller must assure that argv has length at least self.nargs

561 # The bulk of this function is in the (inlined) function process_op.

562 cdef int i, top = -1

563 for i from 0 <= i < self.nops:

564 top = process_op(self.ops[i], self.stack, argv, top)

565 cdef double res = self.stack[0]

566 return res

567

568 def _fast_float_(self, *vars):

569 r"""

570 Returns \code{self} if there are enough arguments, otherwise raises a TypeError.

571

572 EXAMPLES:

573 sage: from sage.ext.fast_eval import fast_float_arg

574 sage: f = fast_float_arg(1)

575 sage: f._fast_float_('x','y') is f

576 True

577 sage: f._fast_float_('x') is f

578 Traceback (most recent call last):

579 ...

580 TypeError: Needs at least 2 arguments (1 provided)

581 """

582 if self.nargs > len(vars):

583 raise TypeError, "Needs at least %s arguments (%s provided)" % (self.nargs, len(vars))

584 return self

585

586 def op_list(self):

587 """

588 Returns a list of string representations of the

589 operations that make up this expression.

590

591 Python and C function calls may be only available by function pointer addresses.

592

593 EXAMPLES:

594 sage: from sage.ext.fast_eval import fast_float_constant, fast_float_arg

595 sage: a = fast_float_constant(17)

596 sage: x = fast_float_arg(0)

597 sage: a.op_list()

598 ['push 17.0']

599 sage: x.op_list()

600 ['load 0']

601 sage: (a*x).op_list()

602 ['push 17.0', 'load 0', 'mul']

603 sage: (a+a*x^2).sqrt().op_list()

604 ['push 17.0', 'push 17.0', 'load 0', 'dup', 'mul', 'mul', 'add', 'call sqrt(1)']

605 """

606 cdef int i

607 return [op_to_string(self.ops[i]) for i from 0 <= i < self.nops]

608

609 def __iter__(self):

610 """

611 Returns the list of operations of self.

612

613 EXAMPLES:

614 sage: from sage.ext.fast_eval import fast_float_arg

615 sage: f = fast_float_arg(0)*2 + 3

616 sage: list(f)

617 ['load 0', 'push 2.0', 'mul', 'push 3.0', 'add']

618 """

619 return iter(self.op_list())

620

621 cpdef bint is_pure_c(self):

622 """

623 Returns True if this function can be evaluated without

624 any python calls (at any level).

625

626 EXAMPLES:

627 sage: from sage.ext.fast_eval import fast_float_constant, fast_float_arg, fast_float_func

628 sage: fast_float_constant(2).is_pure_c()

629 True

630 sage: fast_float_arg(2).sqrt().sin().is_pure_c()

631 True

632 sage: fast_float_func(lambda _: 2).is_pure_c()

633 False

634 """

635 cdef int i

636 for i from 0 <= i < self.nops:

637 if self.ops[i].type == PY_FUNC:

638 return 0

639 return 1

640

641 def python_calls(self):

642 """

643 Returns a list of all python calls used by function.

644

645 EXAMPLES:

646 sage: from sage.ext.fast_eval import fast_float_func, fast_float_arg

647 sage: x = fast_float_arg(0)

648 sage: f = fast_float_func(hash, sqrt(x))

649 sage: f.op_list()

650 ['load 0', 'call sqrt(1)', 'py_call <built-in function hash>(1)']

651 sage: f.python_calls()

652 [<built-in function hash>]

653 """

654 L = []

655 cdef int i

656 for i from 0 <= i < self.nops:

657 if self.ops[i].type == PY_FUNC:

658 L.append((<object>self.ops[i].params.func)[1])

659 return L

660

661 ###################################################################

662 # Basic Arithmatic

663 ###################################################################

664

665 def __add__(left, right):

666 """

667 EXAMPLES:

668 sage: from sage.ext.fast_eval import fast_float_arg

669 sage: f = fast_float_arg(0) + fast_float_arg(1)

670 sage: f(3,4)

671 7.0

672 """

673 return binop(left, right, ADD)

674

675 def __sub__(left, right):

676 """

677 EXAMPLES:

678 sage: from sage.ext.fast_eval import fast_float_arg

679 sage: f = fast_float_arg(0) - fast_float_arg(2)

680 sage: f(3,4,5)

681 -2.0

682 """

683 return binop(left, right, SUB)

684

685 def __mul__(left, right):

686 """

687 EXAMPLES:

688 sage: from sage.ext.fast_eval import fast_float_arg

689 sage: f = fast_float_arg(0) * 2

690 sage: f(17)

691 34.0

692 """

693 return binop(left, right, MUL)

694

695 def __div__(left, right):

696 """

697 EXAMPLES:

698 sage: from sage.ext.fast_eval import fast_float_arg

699 sage: f = fast_float_arg(0) / 7

700 sage: f(14)

701 2.0

702 """

703 return binop(left, right, DIV)

704

705 def __pow__(FastDoubleFunc left, right, dummy):

706 """

707 EXAMPLES:

708 sage: from sage.ext.fast_eval import FastDoubleFunc

709 sage: f = FastDoubleFunc('arg', 0)^2

710 sage: f(2)

711 4.0

712 sage: f = FastDoubleFunc('arg', 0)^4

713 sage: f(2)

714 16.0

715 sage: f = FastDoubleFunc('arg', 0)^-3

716 sage: f(2)

717 0.125

718 sage: f = FastDoubleFunc('arg', 0)^FastDoubleFunc('arg', 1)

719 sage: f(5,3)

720 125.0

721 """

722 if not isinstance(right, FastDoubleFunc):

723 if right == int(right):

724 if right == 1:

725 return left

726 elif right == 2:

727 return left.unop(DUP).unop(MUL)

728 elif right == 3:

729 return left.unop(DUP).unop(DUP).unop(MUL).unop(MUL)

730 elif right == 4:

731 return left.unop(DUP).unop(MUL).unop(DUP).unop(MUL)

732 elif right < 0:

733 return (~left)**(-right)

734 right = FastDoubleFunc('const', right)

735 cdef FastDoubleFunc feval = binop(left, right, TWO_ARG_FUNC)

736 feval.ops[feval.nops-1].params.func = &pow

737 return feval

738

739 def __neg__(FastDoubleFunc self):

740 """

741 EXAMPLE:

742 sage: from sage.ext.fast_eval import fast_float_arg

743 sage: f = -fast_float_arg(0)

744 sage: f(3.5)

745 -3.5

746 """

747 return self.unop(NEG)

748

749 def __abs__(FastDoubleFunc self):

750 """

751 EXAMPLE:

752 sage: from sage.ext.fast_eval import fast_float_arg

753 sage: f = abs(fast_float_arg(0))

754 sage: f(-3)

755 3.0

756 """

757 return self.unop(ABS)

758

759 def abs(FastDoubleFunc self):

760 """

761 EXAMPLE:

762 sage: from sage.ext.fast_eval import fast_float_arg

763 sage: f = fast_float_arg(0).abs()

764 sage: f(3)

765 3.0

766 """

767 return self.unop(ABS)

768

769 def __invert__(FastDoubleFunc self):

770 """

771 EXAMPLE:

772 sage: from sage.ext.fast_eval import fast_float_arg

773 sage: f = ~fast_float_arg(0)

774 sage: f(4)

775 0.25

776 """

777 return self.unop(INVERT)

778

779 def sqrt(self):

780 """

781 EXAMPLE:

782 sage: from sage.ext.fast_eval import fast_float_arg

783 sage: f = fast_float_arg(0).sqrt()

784 sage: f(4)

785 2.0

786 """

787 return self.cfunc(&sqrt)

788

789 ###################################################################

790 # Exponential and log

791 ###################################################################

792

793 def log(self, base=None):

794 """

795 EXAMPLE:

796 sage: from sage.ext.fast_eval import fast_float_arg

797 sage: f = fast_float_arg(0).log()

798 sage: f(2)

799 0.693147180559945...

800 sage: f = fast_float_arg(0).log(2)

801 sage: f(2)

802 1.0

803 sage: f = fast_float_arg(0).log(3)

804 sage: f(9)

805 2.0

806 """

807 if base is None:

808 return self.cfunc(&log)

809 elif base == 2:

810 return self.cfunc(&log2)

811 elif base == 10:

812 return self.cfunc(&log10)

813 else:

814 try:

815 base = fast_float_constant(log(float(base)))

816 except TypeError, e:

817 base = fast_float(base.log())

818 return binop(self.cfunc(&log), base, DIV)

819

820 def exp(self):

821 """

822 EXAMPLE:

823 sage: from sage.ext.fast_eval import fast_float_arg

824 sage: f = fast_float_arg(0).exp()

825 sage: f(1)

826 2.7182818284590451

827 sage: f(100)

828 2.6881171418161356e+43

829 """

830 return self.cfunc(&exp)

831

832 ###################################################################

833 # Rounding

834 ###################################################################

835

836 def ceil(self):

837 """

838 EXAMPLE:

839 sage: from sage.ext.fast_eval import fast_float_arg

840 sage: f = fast_float_arg(0).ceil()

841 sage: f(1.5)

842 2.0

843 sage: f(-1.5)

844 -1.0

845 """

846 return self.cfunc(&ceil)

847

848 def floor(self):

849 """

850 EXAMPLE:

851 sage: from sage.ext.fast_eval import fast_float_arg

852 sage: f = fast_float_arg(0).floor()

853 sage: f(11.5)

854 11.0

855 sage: f(-11.5)

856 -12.0

857 """

858 return self.cfunc(&floor)

859

860 ###################################################################

861 # Trigonometric

862 ###################################################################

863

864 def sin(self):

865 """

866 EXAMPLE:

867 sage: from sage.ext.fast_eval import fast_float_arg

868 sage: f = fast_float_arg(0).sin()

869 sage: f(pi/2)

870 1.0

871 """

872 return self.cfunc(&sin)

873

874 def cos(self):

875 """

876 EXAMPLE:

877 sage: from sage.ext.fast_eval import fast_float_arg

878 sage: f = fast_float_arg(0).cos()

879 sage: f(0)

880 1.0

881 """

882 return self.cfunc(&cos)

883

884 def tan(self):

885 """

886 EXAMPLE:

887 sage: from sage.ext.fast_eval import fast_float_arg

888 sage: f = fast_float_arg(0).tan()

889 sage: f(pi/3)

890 1.73205080756887...

891 """

892 return self.cfunc(&tan)

893

894 def csc(self):

895 """

896 EXAMPLE:

897 sage: from sage.ext.fast_eval import fast_float_arg

898 sage: f = fast_float_arg(0).csc()

899 sage: f(pi/2)

900 1.0

901 """

902 return ~self.sin()

903

904 def sec(self):

905 """

906 EXAMPLE:

907 sage: from sage.ext.fast_eval import fast_float_arg

908 sage: f = fast_float_arg(0).sec()

909 sage: f(pi)

910 -1.0

911 """

912 return ~self.cos()

913

914 def cot(self):

915 """

916 EXAMPLE:

917 sage: from sage.ext.fast_eval import fast_float_arg

918 sage: f = fast_float_arg(0).cot()

919 sage: f(pi/4)

920 1.0...

921 """

922 return ~self.tan()

923

924 def arcsin(self):

925 """

926 EXAMPLE:

927 sage: from sage.ext.fast_eval import fast_float_arg

928 sage: f = fast_float_arg(0).arcsin()

929 sage: f(0.5)

930 0.5235987755982989...

931 """

932 return self.cfunc(&asin)

933

934 def arccos(self):

935 """

936 EXAMPLE:

937 sage: from sage.ext.fast_eval import fast_float_arg

938 sage: f = fast_float_arg(0).arccos()

939 sage: f(sqrt(3)/2)

940 0.5235987755982989...

941 """

942 return self.cfunc(&acos)

943

944 def arctan(self):

945 """

946 EXAMPLE:

947 sage: from sage.ext.fast_eval import fast_float_arg

948 sage: f = fast_float_arg(0).arctan()

949 sage: f(1)

950 0.785398163397448...

951 """

952 return self.cfunc(&atan)

953

954 ###################################################################

955 # Hyperbolic

956 ###################################################################

957

958 def sinh(self):

959 """

960 EXAMPLE:

961 sage: from sage.ext.fast_eval import fast_float_arg

962 sage: f = fast_float_arg(0).sinh()

963 sage: f(log(2))

964 0.75

965 """

966 return self.cfunc(&sinh)

967

968 def cosh(self):

969 """

970 EXAMPLE:

971 sage: from sage.ext.fast_eval import fast_float_arg

972 sage: f = fast_float_arg(0).cosh()

973 sage: f(log(2))

974 1.25

975 """

976 return self.cfunc(&cosh)

977

978 def tanh(self):

979 """

980 EXAMPLE:

981 sage: from sage.ext.fast_eval import fast_float_arg

982 sage: f = fast_float_arg(0).tanh()

983 sage: f(0)

984 0.0

985 """

986 return self.cfunc(&tanh)

987

988 def arcsinh(self):

989 """

990 EXAMPLE:

991 sage: from sage.ext.fast_eval import fast_float_arg

992 sage: f = fast_float_arg(0).arcsinh()

993 sage: f(sinh(5))

994 5.0

995 """

996 return self.cfunc(&asinh)

997

998 def arccosh(self):

999 """

1000 EXAMPLE:

1001 sage: from sage.ext.fast_eval import fast_float_arg

1002 sage: f = fast_float_arg(0).arccosh()

1003 sage: f(cosh(5))

1004 5.0

1005 """

1006 return self.cfunc(&acosh)

1007

1008 def arctanh(self):

1009 """

1010 EXAMPLE:

1011 sage: from sage.ext.fast_eval import fast_float_arg

1012 sage: f = fast_float_arg(0).arctanh()

1013 sage: abs(f(tanh(0.5)) - 0.5) < 0.0000001

1014 True

1015 """

1016 return self.cfunc(&atanh)

1017

1018 cdef FastDoubleFunc cfunc(FastDoubleFunc self, void* func):

1019 cdef FastDoubleFunc feval = self.unop(ONE_ARG_FUNC)

1020 feval.ops[feval.nops - 1].params.func = func

1021 feval.allocate_stack()

1022 return feval

1023

1024 ###################################################################

1025 # Utility functions

1026 ###################################################################

1027

1028 cdef FastDoubleFunc unop(FastDoubleFunc self, char type):

1029 cdef FastDoubleFunc feval = PY_NEW(FastDoubleFunc)

1030 feval.nargs = self.nargs

1031 feval.nops = self.nops + 1

1032 feval.max_height = self.max_height

1033 if type == DUP:

1034 feval.max_height += 1

1035 feval.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op) * feval.nops)

1036 memcpy(feval.ops, self.ops, sizeof(fast_double_op) * self.nops)

1037 feval.ops[feval.nops - 1].type = type

1038 feval.py_funcs = self.py_funcs

1039 feval.allocate_stack()

1040 return feval

1041

1042 cdef FastDoubleFunc binop(_left, _right, char type):

1043 r"""

1044 Returns a function that calculates left and right on the stack, leaving

1045 their results on the top, and then calls operation \code{type}.

1046

1047 EXAMPLES:

1048 sage: from sage.ext.fast_eval import fast_float_arg

1049 sage: f = fast_float_arg(1)

1050 sage: g = fast_float_arg(2) * 11

1051 sage: f.op_list()

1052 ['load 1']

1053 sage: g.op_list()

1054 ['load 2', 'push 11.0', 'mul']

1055 sage: (f+g).op_list()

1056 ['load 1', 'load 2', 'push 11.0', 'mul', 'add']

1057

1058 Correctly calculates the maximum stack heights and number of arguments:

1059 sage: f.max_height

1060 1

1061 sage: g.max_height

1062 2

1063 sage: (f+g).max_height

1064 3

1065 sage: (g+f).max_height

1066 2

1067

1068 sage: f.nargs

1069 2

1070 sage: g.nargs

1071 3

1072 sage: (f+g).nargs

1073 3

1074 """

1075 cdef FastDoubleFunc left, right

1076 try:

1077 left = _left

1078 except TypeError:

1079 left = fast_float(_left)

1080 try:

1081 right = _right

1082 except TypeError:

1083 right = fast_float(_right)

1084 cdef FastDoubleFunc feval = PY_NEW(FastDoubleFunc)

1085 feval.nargs = max(left.nargs, right.nargs)

1086 feval.nops = left.nops + right.nops + 1

1087 feval.max_height = max(left.max_height, right.max_height+1)

1088 feval.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op) * feval.nops)

1089 memcpy(feval.ops, left.ops, sizeof(fast_double_op) * left.nops)

1090 memcpy(feval.ops + left.nops, right.ops, sizeof(fast_double_op) * right.nops)

1091 feval.ops[feval.nops - 1].type = type

1092 if left.py_funcs is None:

1093 feval.py_funcs = right.py_funcs

1094 elif right.py_funcs is None:

1095 feval.py_funcs = left.py_funcs

1096 else:

1097 feval.py_funcs = left.py_funcs + right.py_funcs

1098 feval.allocate_stack()

1099 return feval

1100

1101

1102 def fast_float_constant(x):

1103 """

1104 Return a fast-to-evaluate constant function.

1105

1106 EXAMPLES:

1107 sage: from sage.ext.fast_eval import fast_float_constant

1108 sage: f = fast_float_constant(-2.75)

1109 sage: f()

1110 -2.75

1111

1112 This is all that goes on under the hood:

1113 sage: fast_float_constant(pi).op_list()

1114 ['push 3.14159265359']

1115 """

1116 return FastDoubleFunc('const', x)

1117

1118 def fast_float_arg(n):

1119 """

1120 Return a fast-to-evaluate argument selector.

1121

1122 INPUT:

1123 n -- the (zero-indexed) argument to select

1124

1125 EXAMPLES:

1126 sage: from sage.ext.fast_eval import fast_float_arg

1127 sage: f = fast_float_arg(0)

1128 sage: f(1,2)

1129 1.0

1130 sage: f = fast_float_arg(1)

1131 sage: f(1,2)

1132 2.0

1133

1134 This is all that goes on under the hood:

1135 sage: fast_float_arg(10).op_list()

1136 ['load 10']

1137 """

1138 return FastDoubleFunc('arg', n)

1139

1140 def fast_float_func(f, *args):

1141 """

1142 Returns a wraper around a python function.

1143

1144 INPUT:

1145 f -- a callable python object

1146 args -- a list of FastDoubleFunc inputs

1147

1148 EXAMPLES:

1149 sage: from sage.ext.fast_eval import fast_float_func, fast_float_arg

1150 sage: f = fast_float_arg(0)

1151 sage: g = fast_float_arg(1)

1152 sage: h = fast_float_func(lambda x,y: x-y, f, g)

1153 sage: h(5, 10)

1154 -5.0

1155

1156 This is all that goes on under the hood:

1157 sage: h.op_list() # random memory address

1158 ['load 0', 'load 1', 'py_call <function <lambda> at 0xb62b230>(2)']

1159 """

1160 return FastDoubleFunc('callable', f, *args)

1161

1162

1163 def fast_float(f, *vars):

1164 """

1165 Tries to create a fast float function out of the

1166 input, if possible.

1167

1168 On failure, returns the input unchanged.

1169

1170 INPUT:

1171 f -- an expression

1172 vars -- the names of the arguments

1173

1174 EXAMPLES:

1175 sage: from sage.ext.fast_eval import fast_float

1176 sage: x,y = var('x,y')

1177 sage: f = fast_float(sqrt(x^2+y^2), 'x', 'y')

1178 sage: f(3,4)

1179 5.0

1180

1181 Specifying the argument names is essential, as fast_float objects

1182 only distinguish between arguments by order.

1183 sage: f = fast_float(x-y, 'x','y')

1184 sage: f(1,2)

1185 -1.0

1186 sage: f = fast_float(x-y, 'y','x')

1187 sage: f(1,2)

1188 1.0

1189 """

1190 if isinstance(f, (tuple, list)):

1191 return tuple([fast_float(x, *vars) for x in f])

1192

1193 cdef int i

1194 for i from 0 <= i < len(vars):

1195 if not PY_TYPE_CHECK(vars[i], str):

1196 v = str(vars[i])

1197 # inexact generators display as 1.00..0*x

1198 if '*' in v:

1199 v = v[v.index('*')+1:]

1200 vars = vars[:i] + (v,) + vars[i+1:]

1201

1202 try:

1203 return f._fast_float_(*vars)

1204 except AttributeError:

1205 pass

1206

1207 try:

1208 return FastDoubleFunc('const', float(f))

1209 except TypeError:

1210 pass

1211

1212 try:

1213 from sage.calculus.calculus import SR

1214 return SR(f)._fast_float_(*vars)

1215 except TypeError:

1216 pass

1217

1218 return f

1219

1220 def is_fast_float(x):

1221 return PY_TYPE_CHECK(x, FastDoubleFunc)