Cégep St-Hyacinthe 19 janvier 2012 system:sage
File: /Users/slabbe/Applications/sage-4.7.1/devel/sage/sage/matrix/matrix_integer_dense.pyx
Type: <type ‘builtin_function_or_method’>
Definition: m.determinant(algorithm=’default’, proof=None, stabilize=2)
Docstring:
Return the determinant of this matrix.
INPUT:
algorithm
- 'default' – use ‘pari’ when number of rows less than 50;
otherwise, use ‘padic’
'padic' - uses a p-adic / multimodular algorithm that relies on code in IML and linbox
'linbox' - calls linbox det (you must set proof=False to use this!)
'ntl' - calls NTL’s det function
'pari' - uses PARI
proof - bool or None; if None use proof.linear_algebra(); only relevant for the padic algorithm.
Note
It would be VERY VERY hard for det to fail even with proof=False.
stabilize - if proof is False, require det to be the same for this many CRT primes in a row. Ignored if proof is True.
ALGORITHM: The p-adic algorithm works by first finding a random vector v, then solving A*x = v and taking the denominator d. This gives a divisor of the determinant. Then we compute \det(A)/d using a multimodular algorithm and the Hadamard bound, skipping primes that divide d.
TIMINGS: This is perhaps the fastest implementation of determinants in the world. E.g., for a 500x500 random matrix with 32-bit entries on a core2 duo 2.6Ghz running OS X, Sage takes 4.12 seconds, whereas Magma takes 62.87 seconds (both with proof False). With proof=True on the same problem Sage takes 5.73 seconds. For another example, a 200x200 random matrix with 1-digit entries takes 4.18 seconds in pari, 0.18 in Sage with proof True, 0.11 in Sage with proof False, and 0.21 seconds in Magma with proof True and 0.18 in Magma with proof False.
EXAMPLES:
sage: A = matrix(ZZ,8,8,[3..66]) sage: A.determinant() 0sage: A = random_matrix(ZZ,20,20) sage: D1 = A.determinant() sage: A._clear_cache() sage: D2 = A.determinant(algorithm='ntl') sage: D1 == D2 TrueNext we try the Linbox det. Note that we must have proof=False.
sage: A = matrix(ZZ,5,[1,2,3,4,5,4,6,3,2,1,7,9,7,5,2,1,4,6,7,8,3,2,4,6,7]) sage: A.determinant(algorithm='linbox') Traceback (most recent call last): ... RuntimeError: you must pass the proof=False option to the determinant command to use LinBox's det algorithm sage: A.determinant(algorithm='linbox',proof=False) -21 sage: A._clear_cache() sage: A.determinant() -21A bigger example:
sage: A = random_matrix(ZZ,30) sage: d = A.determinant() sage: A._clear_cache() sage: A.determinant(algorithm='linbox',proof=False) == d True
Latex : $\alpha^2+x^2+\frac{\pi}{6}$
{{{id=190| /// }}}Dessins 2D
{{{id=3| f = sin(1/x) P = plot(f, -10, 10, color='red') P ///by Marshall Hampton.
{{{id=46| %hide npi = RDF(pi) from math import cos,sin def rot(t): return matrix([[cos(t),sin(t)],[-sin(t),cos(t)]]) def pursuit(n,x0,y0,lamb,steps = 100, threshold = .01): paths = [[[x0,y0]]] for i in range(1,n): rx,ry = list(rot(2*npi*i/n)*vector([x0,y0])) paths.append([[rx,ry]]) oldpath = [x[-1] for x in paths] for q in range(steps): diffs = [[oldpath[(j+1)%n][0]-oldpath[j][0],oldpath[(j+1)%n][1]-oldpath[j][1]] for j in range(n)] npath = [[oldpath[j][0]+lamb*diffs[j][0],oldpath[j][1]+lamb*diffs[j][1]] for j in range(n)] for j in range(n): paths[j].append(npath[j]) oldpath = npath return paths html('by William Stein
{{{id=245| /// }}} {{{id=246| /// }}} {{{id=11| %hide import random def ftree(rows, v, i, F): if len(v) > 0: # add a row to g at the ith level. rows.append(v) w = [] for i in range(len(v)): k, _, _ = v[i] if k is None or is_prime(k): w.append((None,None,None)) else: d = random.choice(divisors(k)[1:-1]) w.append((d,k,i)) e = k//d if e == 1: w.append((None,None)) else: w.append((e,k,i)) if len(w) > len(v): ftree(rows, w, i+1, F) def draw_ftree(rows,font): g = Graphics() for i in range(len(rows)): cur = rows[i] for j in range(len(cur)): e, f, k = cur[j] if not e is None: if is_prime(e): c = (1,0,0) else: c = (0,0,.4) g += text(str(e), (j*2-len(cur),-i), fontsize=font, rgbcolor=c) if not k is None and not f is None: g += line([(j*2-len(cur),-i), ((k*2)-len(rows[i-1]),-i+1)], alpha=0.5) return g @interact def factor_tree(n=100, font=(10, (8..20)), redraw=['Redraw']): n = Integer(n) rows = [] v = [(n,None,0)] ftree(rows, v, 0, factor(n)) show(draw_ftree(rows, font), axes=False) /// }}} {{{id=114| /// }}}by William Stein
{{{id=156| @interact def _(N=(100,(2..2000))): html("$\pi(x)$ and $x/(\log(x)-1)$ for $x < %s$"%N) show(plot(prime_pi, 0, N, rgbcolor='red') + plot(x/(log(x)-1), 5, N, rgbcolor='blue')) /// }}}by William Stein
{{{id=115| %hide import urllib class Day: def __init__(self, date, open, high, low, close, volume): self.date = date self.open=float(open); self.high=float(high); self.low=float(low); self.close=float(close) self.volume=int(volume) def __repr__(self): return '%10s %4.2f %4.2f %4.2f %4.2f %10d'%(self.date, self.open, self.high, self.low, self.close, self.volume) class Stock: def __init__(self, symbol): self.symbol = symbol.upper() def __repr__(self): return "%s (%s)"%(self.symbol, self.yahoo()['price']) def yahoo(self): url = 'http://finance.yahoo.com/d/quotes.csv?s=%s&f=%s' % (self.symbol, 'l1c1va2xj1b4j4dyekjm3m4rr5p5p6s7') values = urllib.urlopen(url).read().strip().strip('"').split(',') data = {} data['price'] = values[0] data['change'] = values[1] data['volume'] = values[2] data['avg_daily_volume'] = values[3] data['stock_exchange'] = values[4] data['market_cap'] = values[5] data['book_value'] = values[6] data['ebitda'] = values[7] data['dividend_per_share'] = values[8] data['dividend_yield'] = values[9] data['earnings_per_share'] = values[10] data['52_week_high'] = values[11] data['52_week_low'] = values[12] data['50day_moving_avg'] = values[13] data['200day_moving_avg'] = values[14] data['price_earnings_ratio'] = values[15] data['price_earnings_growth_ratio'] = values[16] data['price_sales_ratio'] = values[17] data['price_book_ratio'] = values[18] data['short_ratio'] = values[19] return data def historical(self): try: return self.__historical except AttributeError: pass symbol = self.symbol def get_data(exchange): name = get_remote_file('http://finance.google.com/finance/historical?q=%s:%s&output=csv'%(exchange, symbol.upper()), verbose=False) return open(name).read() R = get_data('NASDAQ') if "Bad Request" in R: R = get_data("NYSE") R = R.splitlines() headings = R[0].split(',') self.__historical = [] try: for x in reversed(R[1:]): date, opn, high, low, close, volume = x.split(',') self.__historical.append(Day(date, opn,high,low,close,volume)) except ValueError: pass self.__historical = Sequence(self.__historical,cr=True,universe=lambda x:x) return self.__historical def plot_average(self, spline_samples=10): d = self.historical() if len(d) == 0: return text('no historical data at Google Finance about %s'%self.symbol, (0,3)) avg = list(enumerate([(z.high+z.low)/2 for z in d])) P = line(avg) + points(avg, rgbcolor='black', pointsize=4) + \ text(self.symbol, (len(d)*1.05, d[-1].low), horizontal_alignment='right', rgbcolor='black') if spline_samples > 0: k = 250//spline_samples spl = spline([avg[i*k] for i in range(len(d)//k)] + [avg[-1]]) P += plot(spl, (0,len(d)+30), color=(0.7,0.7,0.7)) P.xmax(260) return P def plot_diff(self): d = self.historical() if len(d) == 0: return text('no historical data at Google Finance about %s'%self.symbol, (0,3)) diff = [] for i in range(1, len(d)): z1 = d[i]; z0 = d[i-1] diff.append((i, (z1.high+z1.low)/2 - (z0.high + z0.low)/2)) P = line(diff,thickness=0.5) + points(diff, rgbcolor='black', pointsize=4) + \ text(self.symbol, (len(d)*1.05, 0), horizontal_alignment='right', rgbcolor='black') P.xmax(260) return P symbols = ['bsc', 'vmw', 'sbux', 'aapl', 'amzn', 'goog', 'wfmi', 'msft', 'yhoo', 'ebay', 'java', 'rht', ]; symbols.sort() stocks = dict([(s,Stock(s)) for s in symbols]) @interact def data(symbol = symbols, other_symbol='', spline_samples=(8,[0..15])): if other_symbol != '': symbol = other_symbol S = Stock(symbol) html('%s | %s |
Animations
{{{id=20| reset() a = animate([sin(x + float(k)) for k in srange(0,2*3.14,0.31)], xmin=0, xmax=2*pi, figsize=[2,1]) /// }}} {{{id=19| a.show() # nécessite la commande convert /// }}} {{{id=187| /// }}}Quantumino Solver
Bientôt, on pourra résoudre le Quantumino Puzzle avec Sage. Voir ce vidéo sur youtube. Le ticket #11379 en question a été fusionné à Sage-4.7.2.alpha2 en août 2011.
{{{id=186| from sage.games.quantumino import show_pentaminos show_pentaminos() /// }}} {{{id=68| from sage.games.quantumino import QuantuminoSolver s = QuantuminoSolver(5).solve().next() s.show3d() /// }}} {{{id=244| s /// Quantumino state where the following pentamino is put aside : Polyomino: [(0, 0, 0), (1, 0, 0), (1, 0, 1), (1, 1, 0), (2, 0, 1)], Color: red }}} {{{id=185| s.show3d() /// }}}The Sage notebook allows transparently editing and compiling Cython code simply by typing %cython at the top of a cell and evaluate it. Variables and functions defined in a Cython cell are imported into the running session.
Example 1, pure Python
Here is some simple Python code to numerically integrate the function $f(x) = x^2$.
{{{id=47| /// }}} {{{id=1| from math import sin def f(x): return sin(x**2) def integrate_f_py(a, b, N): s = 0 dx = (b-a)/N for i in range(N): s += f(a+i*dx) return s * dx /// }}} {{{id=62| timeit('integrate_f_py(0, 1, 1000)', number=50) /// }}} {{{id=7| /// }}} {{{id=63| /// }}}Example 1, compiled with Cython (no other changes)
Simply compiling this in Cython gives a speedup.
{{{id=5| %cython from math import sin def f(x): return sin(x**2) def integrate_f_cy0(a, b, N): s = 0 dx = (b-a)/N for i in range(N): s += f(a+i*dx) return s * dx /// }}} {{{id=4| timeit('integrate_f_cy0(0, 1, 1000)', number=50) /// }}} {{{id=8| /// }}} {{{id=69| /// }}}Example 1, typed and compiled with Cython
Adding some static type declarations makes a much greater difference.
{{{id=10| %cython from math import sin def f(double x): return sin(x**2) def integrate_f_cy(double a, double b, int N): cdef int i cdef double s, dx s = 0 dx = (b-a)/N for i in range(N): s += f(a+i*dx) return s * dx /// }}} {{{id=72| timeit('integrate_f_cy(0, 1, 1000)') /// }}} {{{id=12| 18500 /489.0 /// }}} {{{id=17| /// }}}Example 2, pure Python
Here is a Python function that computes the sum of the first $n$ positive integers.
{{{id=29| def mysum_py(n): s = 0 for k in range(n): s += k return s /// }}} {{{id=31| time mysum_py(10^6) /// 499999500000 Time: CPU 2.26 s, Wall: 2.27 s }}} {{{id=32| /// }}} {{{id=74| /// }}}Example 2, just compiled with Cython
Simply compiling this function with Cython provides a speedup.
{{{id=83| %cython def mysum_cy0(n): s = 0 for k in range(n): s += k return s /// }}} {{{id=84| time mysum_cy0(10^6) /// 499999500000L Time: CPU 0.24 s, Wall: 0.24 s }}} {{{id=85| 2.09 / 0.22 /// 9.50000000000000 }}} {{{id=86| /// }}}Example 2, typed and compiled with Cython
Adding some static type declarations makes a much greater difference.
{{{id=88| %cython def mysum_cy1(n): cdef int k cdef long long s s = 0 for k in range(n): s += k return s /// }}} {{{id=42| time mysum_cy1(10^6) /// 499999500000L Time: CPU 0.00 s, Wall: 0.00 s }}} {{{id=89| 2.09 / 0.00 /// +infinity }}} {{{id=179| timeit('mysum_cy1(10^6)') /// 625 loops, best of 3: 1.58 ms per loop }}} {{{id=90| 2.09/0.00158 /// 1322.78481012658 }}}Two important (but minor) differences between Sage language and Python
Integer division in Python :
in Sage:
{{{id=125| 2/3 + 4/5 + 1/7 /// }}}Exponent (^) in Python :
{{{id=130| %python 10^14 #exclusive OR /// }}}in Sage :
{{{id=149| 10^14 /// }}}The preparser
{{{id=154| preparse('2/3 + 2^3 + 3.0') /// }}} {{{id=35| preparse('2^3') /// }}} {{{id=163| /// }}}