Upgrade to Pro — share decks privately, control downloads, hide ads and more …

First try for CAS, SymPy with codegen

dv
June 05, 2016

First try for CAS, SymPy with codegen

dv

June 05, 2016
Tweet

More Decks by dv

Other Decks in Programming

Transcript

  1. Who am I ? • 我朋友都叫我 "dv" • Pythonista •

    FOSS 熱愛 • 快畢業的大學生 • 即將進入四個月國軍 Online • 出來後準備找工作
  2. CAS (Computer Algebra System) 一些知名的軟體: • Maxima (GPL) • SageMath

    (GPL) • SymPy (BSD License) • Mathematica (Proprietary) • Magma (Proprietary) • Maple (Proprietary) • MATLAB - Symbolic Math Toolbox (Proprietary)
  3. 什麼是 SymPy? • Python 的 Library • Open Source (BSD

    License) • 針對數學符號做處理 (Symbolic Computation)
  4. 什麼是 SymPy? • Python 的 Library • Open Source (BSD

    License) • 針對數學符號做處理 (Symbolic Computation) • 安裝方便 ( 依賴的軟體相對少 )
  5. Try Jupyter! try.jupyter.org 1. 選 "New" 然後 "Python 3" 2.

    輸入 "import sympy" 然後按 "Shift + Enter" 3. 現在你可以使用任何 SymPy 的功能了
  6. x, y, z = symbols ( "x, y, z" )

    from sympy import symbols
  7. x, y = symbols ( "x, y" ) expr =

    x**2 + 2*x*y + y**2
  8. expr = sympify ( "x**2 + 2*x*y + y**2" )

    from sympy import sympify
  9. expr = sympify ( "x**2 + 2*x*y + y**2" )

    result = expr.subs ( {"x": 1, "y": 2} )
  10. expr = sympify ( "x**2 + 2*x*y + y**2" )

    result = expr.evalf ( subs = {"x": 1, "y": 2} )
  11. expr = sympify ( "x**4 - 1" ) result =

    solve ( expr ) [-1, 1, -I, I] from sympy import solve
  12. expr1 = sympify ( "x + y - 10" )

    expr2 = sympify ( "x - y + 6" ) result = solve ( [expr1, expr2] ) {y: 8, x: 2}
  13. static const double half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */

    S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ FreeBSD - msun (k_sin.c) 方程式係數
  14. GNU C Library (s_sin.c) static const double sn3 = -1.66666666666664880952546298448555E-01,

    sn5 = 8.33333214285722277379541354343671E-03, cs2 = 4.99999999999999999999950396842453E-01, cs4 = -4.16666666666664434524222570944589E-02, cs6 = 1.38888874007937613028114285595617E-03; 方程式係數
  15. musl (__sin.c) static const double S1 = -1.66666666666666324348e-01, /* 0xBFC55555,

    0x55555549 */ S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ 方程式係數
  16. # NumPy universal function ufuncify ( args, expr, backend="numpy" )

    from sympy.utilities.autowrap import ufuncify
  17. # Python & NumPy's C API # and C compiler

    ufuncify ( x, x**x ) <ufunc 'wrapper_module_0'>
  18. # broadcast and type theano_function ( [x], [x**x], dims={x: 1},

    dtype={x: 'float64'}) <theano.compile.function_module.Function at 0x7fd9a0bede48>
  19. poly, err = chebyfit ( sin, [0, pi/4], 10, error=True)

    from mpmath import chebyfit, sin, pi poly: [2.53705546510064e-6, 4.83882390985545e-7, -0.000198913241162381, 2.96736201694264e-7, 0.00833322702020108, 2.28887435232138e-8, -0.166666669474202, 1.73888384312788e-10, 0.999999999995841, 1.63478544118895e-14]
  20. poly, err = chebyfit ( sin, [0, pi/4], 10, error=True)

    from mpmath import chebyfit, sin, pi err: 1.94337098561508e-14
  21. Python 直譯器實作 • CPython ( 官方 ) • PyPy -使用

    Python 的子集來實作,採用 JIT 與不同的 Garbage Collection 技術 • Pyston -基於 LLVM 實作 JIT • Jython -把 Python code 轉成 Java bytecode 跑在 JVM 上 • IronPython -用 .NET Framework 實作 • …
  22. PyPy pypy.org • 專案狀態 ◦ Python 2.7 有支援 ◦ Python

    3 目前指支援到 Python 3.3 (alpha) • 效能比較可以看 speed.pypy.org
  23. Pyston github.com/dropbox/pyston • 由 Dropbox 於 2013 發起的專案 • 使用

    C++11 撰寫 • 基於 LLVM 實作 JIT • 目前方向著重在 Python 2.7
  24. C

  25. codegen ( ("func", x**x), "c", header=False) #include "func.h" #include <math.h>

    double func(double x) { double func_result; func_result = pow(x, x); return func_result; }
  26. SymPy Expression Tree 把 SymPy 的 Expression Tree 轉成 DOT

    語言描述 (接著可以用 Graphviz 這類工具轉成圖檔) from sympy.printing.dot import dotprint
  27. _print ( expr ) _print_Rational ( expr ) _print_Number (

    expr ) _print_AtomicExpr ( expr ) _print_Atom ( expr ) _print_Basic ( expr )
  28. SymPy 內其他 Modules • Calculus • Physics • Lie Algebra

    • Linear Algebra • Number Theory • Differential Geometry • Set Theory • ODE • PDE • Statistics • Logic • …
  29. 基於 SymPy 的專案 • PyDy - Multibody Dynamics with Python

    • SageMath - 有把 SymPy 整合進去 • SfePy - Simple Finite Elements in Python • …