CoCalc Shared Files2017-04-11-112324.sagews
Author: Craig Larson
Views : 6
3+4

7
3==4

False
3+4

7
is_prime(47)

True
is_prime(1111111111111111111)

True
is_prime(12345678910987654321)

True
diff(x^x)

x^x*(log(x) + 1)
for x in [4,5,6,7]:
print x^2

16 25 36 49
pete = graphs.PetersenGraph()

show(pete)

d3-based renderer not yet implemented
pete.order()

10
pete.size()

15
cory = Graph(7)

show(cory)

d3-based renderer not yet implemented
cory.add_edge(0,5)

show(cory)

d3-based renderer not yet implemented
load("conjecturing.py")

objects = [2,3,5,60]

factor(60)

2^2 * 3 * 5
len(factor(60))

3
%lisp

(* (+ 3 4) (- 5 3))

14
def distinct_primes(n):
return len(factor(n))

distinct_primes(60)

3
47.digits()

[7, 4]
sum(47.digits())

11
def sum_digits(n):
return sum(n.digits())

sum_digits(47)

11
def number(n):
return n

number(135)

135
L=[5,6,7]
L.index(6)

1
load("conjecturing.py")

n(ln(10))

2.30258509299405
log(10,10)

1
log(100,10)

2
n(log(10))

2.30258509299405
load("conjecturing.py")

def sum_digits(n):
return sum(n.digits())
def number(n):
return n
def distinct_primes(n):
return len(factor(n))

objects = [2,3,5,60]
invariants = [number, sum_digits, distinct_primes]

conjecture(objects,invariants,invariants.index(number))


[number(x) <= 10^(log(sum_digits(x))/log(10) + 1), number(x) <= sum_digits(x)^distinct_primes(x)]
def conj1(n):
return numerical_approx(10^(log(sum_digits(n))/log(10) + 1))

def conj2(n):
return numerical_approx(sum_digits(n)^distinct_primes(n))

conj2(2)

2.00000000000000
conj2(3)

3.00000000000000
conj2(5)

5.00000000000000
conj2(60)

216.000000000000
n

<function numerical_approx at 0x7feaadc94848>
conj1(2)

20.0000000000000
conj1(3)

30.0000000000000
conj1(5)

50.0000000000000
conj1(60)

60.0000000000000
def number(n):
return n


#47 is a counterexample to conjecture 2
#add to the program and run new conjectures

def sum_digits(n):
return sum(n.digits())
def number(n):
return n
def distinct_primes(n):
return len(factor(n))

objects = [2,3,5,60,47]
invariants = [number, sum_digits, distinct_primes]

conjecture(objects,invariants,invariants.index(number))

[number(x) <= log(log(10^e^sum_digits(x))/log(10)), number(x) <= ceil((sum_digits(x) + 1)^cosh(distinct_primes(x))), number(x) <= 10^(log(sum_digits(x))/log(10) + 1), number(x) <= floor(e^(sum_digits(x) - 1)), number(x) <= floor(cosh(distinct_primes(x))^sum_digits(x))]
3+4

7
x=7
def test(y):
y = y+1
return y

test(x)
x

8 7
L=[1,2,3]
def test2(X):
X[0]=47
return X

test2(L)
L

[47, 2, 3] [47, 2, 3]
sage: import sage.libs.mpmath.all as mpmath
sage: mpmath.mp.dps = 30
sage: mpmath.pslq([sqrt(n) for n in range(2, 8+1)])

[2, 0, 0, 0, 0, 0, -1]
integral(x^2*log(x)/((x^2-1)*(x^4+1)),x,0,1)

1/32*sqrt(2)*pi^2*(sqrt(2) - 1)
numerical_integral(x^2*log(x)/((x^2-1)*(x^4+1)),0,1)

(0.18067126259059388, 8.044808609022901e-11)
mpmath.mp.identify(0.881373587019543)

mpmath.mp.identify(0.22222222222222)

'((36-sqrt(0))/162)'
mpmath.mp.identify?

File: /projects/sage/sage-7.5/local/lib/python2.7/site-packages/mpmath/identification.py
Signature : mpmath.mp.identify(ctx, x, constants=[], tol=None, maxcoeff=1000, full=False, verbose=False)
Docstring :
Given a real number x, "identify(x)" attempts to find an exact
formula for x. This formula is returned as a string. If no match is
found, "None" is returned. With "full=True", a list of matching
formulas is returned.

As a simple example, "identify()" will find an algebraic formula
for the golden ratio:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> identify(phi)
'((1+sqrt(5))/2)'

"identify()" can identify simple algebraic numbers and simple
combinations of given base constants, as well as certain basic
transformations thereof. More specifically, "identify()" looks for
the following:

1. Fractions

2. Quadratic algebraic numbers

3. Rational linear combinations of the base constants

4. Any of the above after first transforming x into f(x)
where f(x) is 1/x, sqrt x, x^2, log x or exp x, either
directly or with x or f(x) multiplied or divided by one of
the base constants

5. Products of fractional powers of the base constants and
small integers

Base constants can be given as a list of strings representing
mpmath expressions ("identify()" will "eval" the strings to
numerical values and use the original strings for the output), or
as a dict of formula:value pairs.

In order not to produce spurious results, "identify()" should be
used with high precision; preferably 50 digits or more.

**Examples**

Simple identifications can be performed safely at standard
precision. Here the default recognition of rational, algebraic, and
exp/log of algebraic numbers is demonstrated:

>>> mp.dps = 15
>>> identify(0.22222222222222222)
'(2/9)'
>>> identify(1.9662210973805663)
'sqrt(((24+sqrt(48))/8))'
>>> identify(4.1132503787829275)
'exp((sqrt(8)/2))'
>>> identify(0.881373587019543)
'log(((2+sqrt(8))/2))'

By default, "identify()" does not recognize pi. At standard
precision it finds a not too useful approximation. At slightly
increased precision, this approximation is no longer accurate
enough and "identify()" more correctly returns "None":

>>> identify(pi)
'(2**(176/117)*3**(20/117)*5**(35/39))/(7**(92/117))'
>>> mp.dps = 30
>>> identify(pi)
>>>

Numbers such as pi, and simple combinations of user-defined
constants, can be identified if they are provided explicitly:

>>> identify(3*pi-2*e, ['pi', 'e'])
'(3*pi + (-2)*e)'

Here is an example using a dict of constants. Note that the
constants need not be "atomic"; "identify()" can just as well
express the given number in terms of expressions given by formulas:

>>> identify(pi+e, {'a':pi+2, 'b':2*e})
'((-2) + 1*a + (1/2)*b)'

Next, we attempt some identifications with a set of base constants.
It is necessary to increase the precision a bit.

>>> mp.dps = 50
>>> base = ['sqrt(2)','pi','log(2)']
>>> identify(0.25, base)
'(1/4)'
>>> identify(3*pi + 2*sqrt(2) + 5*log(2)/7, base)
'(2*sqrt(2) + 3*pi + (5/7)*log(2))'
>>> identify(exp(pi+2), base)
'exp((2 + 1*pi))'
>>> identify(1/(3+sqrt(2)), base)
'((3/7) + (-1/7)*sqrt(2))'
>>> identify(sqrt(2)/(3*pi+4), base)
'sqrt(2)/(4 + 3*pi)'
>>> identify(5**(mpf(1)/3)*pi*log(2)**2, base)
'5**(1/3)*pi*log(2)**2'

An example of an erroneous solution being found when too low
precision is used:

>>> mp.dps = 15
>>> identify(1/(3*pi-4*e+sqrt(8)), ['pi', 'e', 'sqrt(2)'])
'((11/25) + (-158/75)*pi + (76/75)*e + (44/15)*sqrt(2))'
>>> mp.dps = 50
>>> identify(1/(3*pi-4*e+sqrt(8)), ['pi', 'e', 'sqrt(2)'])
'1/(3*pi + (-4)*e + 2*sqrt(2))'

**Finding approximate solutions**

The tolerance "tol" defaults to 3/4 of the working precision.
Lowering the tolerance is useful for finding approximate matches.
We can for example try to generate approximations for pi:

>>> mp.dps = 15
>>> identify(pi, tol=1e-2)
'(22/7)'
>>> identify(pi, tol=1e-3)
'(355/113)'
>>> identify(pi, tol=1e-10)
'(5**(339/269))/(2**(64/269)*3**(13/269)*7**(92/269))'

With "full=True", and by supplying a few base constants, "identify"
can generate almost endless lists of approximations for any number
(the output below has been truncated to show only the first few):

>>> for p in identify(pi, ['e', 'catalan'], tol=1e-5, full=True):
...     print(p)
...  # doctest: +ELLIPSIS
e/log((6 + (-4/3)*e))
(3**3*5*e*catalan**2)/(2*7**2)
sqrt(((-13) + 1*e + 22*catalan))
log(((-6) + 24*e + 4*catalan)/e)
exp(catalan*((-1/5) + (8/15)*e))
catalan*(6 + (-6)*e + 15*catalan)
sqrt((5 + 26*e + (-3)*catalan))/e
e*sqrt(((-27) + 2*e + 25*catalan))
log(((-1) + (-11)*e + 59*catalan))
((3/20) + (21/20)*e + (3/20)*catalan)
...

The numerical values are roughly as close to pi as permitted by
the specified tolerance:

>>> e/log(6-4*e/3)
3.14157719846001
>>> 135*e*catalan**2/98
3.14166950419369
>>> sqrt(e-13+22*catalan)
3.14158000062992
>>> log(24*e-6+4*catalan)-1
3.14158791577159

**Symbolic processing**

The output formula can be evaluated as a Python expression. Note
however that if fractions (like '2/3') are present in the formula,
Python's "eval()" may erroneously perform integer division. Note
also that the output is not necessarily in the algebraically
simplest form:

>>> identify(sqrt(2))
'(sqrt(8)/2)'

As a solution to both problems, consider using SymPy's "sympify()"
to convert the formula into a symbolic expression. SymPy can be
used to pretty-print or further simplify the formula symbolically:

>>> from sympy import sympify # doctest: +SKIP
>>> sympify(identify(sqrt(2))) # doctest: +SKIP
2**(1/2)

Sometimes "identify()" can simplify an expression further than a
symbolic algorithm:

>>> from sympy import simplify # doctest: +SKIP
>>> x = sympify('-1/(-3/2+(1/2)*5**(1/2))*(3/2-1/2*5**(1/2))**(1/2)') # doctest: +SKIP
>>> x # doctest: +SKIP
(3/2 - 5**(1/2)/2)**(-1/2)
>>> x = simplify(x) # doctest: +SKIP
>>> x # doctest: +SKIP
2/(6 - 2*5**(1/2))**(1/2)
>>> mp.dps = 30 # doctest: +SKIP
>>> x = sympify(identify(x.evalf(30))) # doctest: +SKIP
>>> x # doctest: +SKIP
1/2 + 5**(1/2)/2

(In fact, this functionality is available directly in SymPy as the
function "nsimplify()", which is essentially a wrapper for
"identify()".)

**Miscellaneous issues and limitations**

The input x must be a real number. All base constants must be
positive real numbers and must not be rationals or rational linear
combinations of each other.

The worst-case computation time grows quickly with the number of
base constants. Already with 3 or 4 base constants, "identify()"
may require several seconds to finish. To search for relations
among a large number of constants, you should consider using
"pslq()" directly.

The extended transformations are applied to x, not the constants
separately. As a result, "identify" will for example be able to
recognize "exp(2*pi+3)" with "pi" given as a base constant, but not
"2*exp(pi)+3". It will be able to recognize the latter if "exp(pi)"
is given explicitly as a base constant.


#a function that returns the first digit of an integer
def first_digit(n):
return n.digits().pop()

#a function that counts the first digitso of each kind from a list of numbers L
def first_digit_counts(L):
data = [0]*10
for x in L:
#print x, first_digit(x), data[first_digit(x)]
first = first_digit(x)
data[first] = data[first] + 1
return data

#a function that gives the (empirical) probabilities that the integer from a list L starts with digit d
def first_digit_distribution(L):
counts = first_digit_counts(L)
length = len(L)
return [n(counts[i]/length, digits=3) for i in [0..9]]

counts=first_digit_counts(srange(1,1000000000))

bar_chart(counts)

mpmath.mp.identify(0.180671262590654)

3+4

7
@interact
def _(f=x^2, a=-3, b=3):
show(plot(f,(x,a,b)))

Interact: please open in CoCalc
@interact
def _(f=input_box(x^2,width=20),
color=color_selector(widget='colorpicker', label=""),
axes=True,
fill=True,
zoom=range_slider(-3,3,default=(-3,3))):
show(plot(f,(x,zoom[0], zoom[1]), color=color, axes=axes,fill=fill))

@interact
def _(f=input_box(x^2,width=20),
color=color_selector(widget='colorpicker', label=""),
axes=True,
fill=True,
zoom=range_slider(-3,3,default=(-3,3))):
show(plot(f,(x,zoom[0], zoom[1]), color=color, axes=axes,fill=fill))

3+4

7
@interact
def _(f=input_box(x^2,width=20),
axes=True,
fill=True,
zoom=range_slider(-3,3,default=(-3,3))):
show(plot(f,(x,zoom[0], zoom[1]), axes=axes,fill=fill))

@interact
def _(f=input_box(x^2,width=20),
axes=True,
fill=True,
zoom=range_slider(-3,3,default=(-3,3))):
show(plot(f,(x,zoom[0], zoom[1]), axes=axes,fill=fill))

Interact: please open in CoCalc
@interact
def _(m=('matrix', identity_matrix(2)), auto_update=False):
print m.eigenvalues()

Interact: please open in CoCalc
var("z")
complex_plot(z**2, (-5, 5), (-5, 5))

z
complex_plot?

File: /projects/sage/sage-7.5/src/sage/misc/lazy_import.pyx
Signature : complex_plot(f, xrange, yrange, plot_points=100, interpolation='catrom', **kwds)
Docstring :
"complex_plot" takes a complex function of one variable, f(z) and
plots output of the function over the specified "xrange" and
"yrange" as demonstrated below. The magnitude of the output is
indicated by the brightness (with zero being black and infinity
being white) while the argument is represented by the hue (with red
being positive real, and increasing through orange, yellow, ... as
the argument increases).

"complex_plot(f, (xmin, xmax), (ymin, ymax), ...)"

INPUT:

* "f" -- a function of a single complex value x + iy

* "(xmin, xmax)" -- 2-tuple, the range of "x" values

* "(ymin, ymax)" -- 2-tuple, the range of "y" values

The following inputs must all be passed in as named parameters:

* "plot_points" -- integer (default: 100); number of points to
plot in each direction of the grid

* "interpolation" -- string (default: "'catrom'"), the
interpolation method to use: "'bilinear'", "'bicubic'",
"'spline16'", "'spline36'", "'quadric'", "'gaussian'", "'sinc'",
"'bessel'", "'mitchell'", "'lanczos'", "'catrom'", "'hermite'",
"'hanning'", "'hamming'", "'kaiser'"

EXAMPLES:

Here we plot a couple of simple functions:

sage: complex_plot(sqrt(x), (-5, 5), (-5, 5))
Graphics object consisting of 1 graphics primitive

sage: complex_plot(sin(x), (-5, 5), (-5, 5))
Graphics object consisting of 1 graphics primitive

sage: complex_plot(log(x), (-10, 10), (-10, 10))
Graphics object consisting of 1 graphics primitive

sage: complex_plot(exp(x), (-10, 10), (-10, 10))
Graphics object consisting of 1 graphics primitive

A function with some nice zeros and a pole:

sage: f(z) = z^5 + z - 1 + 1/z
sage: complex_plot(f, (-3, 3), (-3, 3))
Graphics object consisting of 1 graphics primitive

Here is the identity, useful for seeing what values map to what
colors:

sage: complex_plot(lambda z: z, (-3, 3), (-3, 3))
Graphics object consisting of 1 graphics primitive

The Riemann Zeta function:

sage: complex_plot(zeta, (-30,30), (-30,30))
Graphics object consisting of 1 graphics primitive

Extra options will get passed on to show(), as long as they are
valid:

sage: complex_plot(lambda z: z, (-3, 3), (-3, 3), figsize=[1,1])
Graphics object consisting of 1 graphics primitive

sage: complex_plot(lambda z: z, (-3, 3), (-3, 3)).show(figsize=[1,1]) # These are equivalent

@interact
def julia_plot(expo = slider(-10,10,0.1,2), \
iterations=slider(1,100,1,30), \
c_real = slider(-2,2,0.01,0.5), \
c_imag = slider(-2,2,0.01,0.5), \
zoom_x = range_slider(-2,2,0.01,(-1.5,1.5)), \
zoom_y = range_slider(-2,2,0.01,(-1.5,1.5))):
var('z')
I = CDF.gen()
f(z) = z^expo + c_real + c_imag*I
ff_j = fast_callable(f, vars=[z], domain=CDF)

def julia(z):
for i in range(iterations):
z = ff_j(z)
if abs(z) > 2:
return z
return z
print 'z <- z^%s + (%s+%s*I)' % (expo, c_real, c_imag)
complex_plot(julia, zoom_x,zoom_y, plot_points=200, dpi=150).show(frame=True, aspect_ratio=1)

Interact: please open in CoCalc

@interact
def _(n=(5..100)):
Poset(([1..n], lambda x, y: y%x == 0) ).show()

Interact: please open in CoCalc
5//2

2
matrix_plot(matrix(2,[1,2,3,4]))

def sierpinski(N):
'''Generates the Sierpinski triangle by taking the modulo-2
of each element in Pascal's triangle'''
S=[([0]*(N//2-a//2))+
[binomial(a,b)%2 for b in range(a+1)]+
([0]*(N//2-a//2)) for a in range(0,N,2)]
return S

@interact
def _(N=slider([2 ** a for a in range(12)],
label='Number of iterations', default=64),
size=slider(1, 20, label='Size', step_size=1, default=9)):
M = sierpinski(2 * N)
matrix_plot(M, cmap='binary').show(figsize=[size, size])

Interact: please open in CoCalc
@interact
def rwalk3d(n=slider(50,1000,step_size=1), frame=True):
pnt = [0,0,0]
v = [copy(pnt)]
for i in range(n):
pnt[0] += random()-0.5
pnt[1] += random()-0.5
pnt[2] += random()-0.5
v.append(copy(pnt))
show(line3d(v,color='black'),aspect_ratio=[1,1,1],frame=frame)

Interact: please open in CoCalc
e

e e
n(e)

2.71828182845905