rationalize

# By Ryan Ruff (c) Aug 12, 2010 # http://twitter.com/suburbanlion # http://suburbanlion.com # Released under the terms of the CC-BY License # http://creativecommons.org/licenses/by/3.0/ # You are free to use, share and remix this work, provided the source is attributed and # all modifications are released under a similar license # This program is meant to solve a problem posed by G. Michael Guy (@GMichaelGuy) via Twitter # http://twitter.com/GMichaelGuy/status/20763608488 # Suppose we want to rationalize the expression given by: 1/(7+5^(1/3)+3^(2/3)+2^(1/3)) # One way to look at this problem is to think of it as the dot product of some vector 'x' # with some basis vector 'b' generated from 2^(1/3),3^(1/3),5^(1/3) # This framework allows use to solve for the rationalization by using linear algebra # import operators for addition and multiplication from operator import add from operator import mul # specify max root to use (if varied, using the LCM should work) nthRoot = 3 # define the primes used for our basis vector p = [2,3,5] # count how many primes we have numPrimes = 0 for i in p: numPrimes = numPrimes + 1 # helper function to serialize the numPrimes*nthRoot possible combinations of powers def powers(n): o = range(numPrimes); t = n; for i in range(numPrimes): o[i] = t % nthRoot; t = (t - (t % nthRoot))/nthRoot; return o; # store this vector so we only compute it once powerVector = map(powers,range(numPrimes^nthRoot)) # generate a basis vector from our serialization b = vector([reduce(mul,map(lambda j,k:j^(k/nthRoot),p,i)) for i in powerVector]) print("Basis Vector (b):") print(b) # define the denominator of our expression by taking the coefficients of that expression relative to basis 'b' # to produce our vector x x = vector([7,1,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0]) print("Expression vector (x):") print(x) # verify that dot(x,B) gives us the expression in problem print("<x,b>:") print(x*b) # we want to find some vector y such that <x,b>*<y,b> = 1 # in order to do that, we want to first find some matrix M such that b*M*y = <x,b>*<y,b> # this is where it gets a little tricky # define a helper function to compute the sum of two vectors modulo 'nthRoot' def add_mod_n(u,v): return map(lambda a, b: (a+b)%nthRoot,u,v); # define a helper function to compute the rational part of the product of two basis vectors def mult_coef(c,r): modSum = add_mod_n(powerVector[c],powerVector[r]); o = 1; for i in range(numPrimes): o = o*p[i]^((powerVector[c][i]+powerVector[r][i]-modSum[i])/3); return o; # define a helper function for finding the conjugate of some basis vector # i.e. the product of return value and c-th component of basis vector lies along r-th component of basis vector def get_conj(c,r): u = map(lambda i,j: (nthRoot+i-j)%nthRoot,powerVector[r],powerVector[c]); return reduce(add, map(lambda i,j: (3^j)*i,u,range(numPrimes))); # generate a row of our product matrix def mult_row(n, v, r): return [v[get_conj(c,r) ]*mult_coef(c,get_conj(c,r)) for c in range(n)]; # generate matrix that simulates multiplication by <b,v> def mult_mat(v): return [mult_row(numPrimes^nthRoot,v,i) for i in range(numPrimes^nthRoot)]; # create our matrix for multiplying by <b,x> M = matrix(mult_mat(x)) # now all thats left is to invert that matrix and apply it to a oneVector oneVec = vector([1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0]) y = M \ oneVec print("y:") print(y) # out put this in a more conventional form by dotting with basis print("<b,y>:") print(b*y) # compare this result with Wolfram Alpha # http://www.wolframalpha.com/input/?i=1/(7%2B5^(1/3)%2B3^(2/3)%2B2^(1/3)) 

Basis Vector (b):
(1, 2^(1/3), 2^(2/3), 3^(1/3), 2^(1/3)*3^(1/3), 2^(2/3)*3^(1/3),
3^(2/3), 2^(1/3)*3^(2/3), 2^(2/3)*3^(2/3), 5^(1/3), 2^(1/3)*5^(1/3),
2^(2/3)*5^(1/3), 3^(1/3)*5^(1/3), 2^(1/3)*3^(1/3)*5^(1/3),
2^(2/3)*3^(1/3)*5^(1/3), 3^(2/3)*5^(1/3), 2^(1/3)*3^(2/3)*5^(1/3),
2^(2/3)*3^(2/3)*5^(1/3), 5^(2/3), 2^(1/3)*5^(2/3), 2^(2/3)*5^(2/3),
3^(1/3)*5^(2/3), 2^(1/3)*3^(1/3)*5^(2/3), 2^(2/3)*3^(1/3)*5^(2/3),
3^(2/3)*5^(2/3), 2^(1/3)*3^(2/3)*5^(2/3), 2^(2/3)*3^(2/3)*5^(2/3))
Expression vector (x):
(7, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0)
<x,b>:
2^(1/3) + 3^(2/3) + 5^(1/3) + 7
y:
(6085751511148292591/44253206255743824031,
-113737199574721883/6321886607963403433,
14159343943246004/6321886607963403433,
45903962039615943/6321886607963403433,
-119414810568569283/44253206255743824031,
7609214068173315/12643773215926806866,
-118070567437573862/6321886607963403433,
30694440328632370/6321886607963403433,
-79218868344463653/88506412511487648062,
-114594372661917362/6321886607963403433,
27999438968454556/6321886607963403433,
-63303351032438961/88506412511487648062,
-135774633227300964/44253206255743824031,
10778499819152451/6321886607963403433,
-7004738263654053/12643773215926806866,
32291803150446358/6321886607963403433,
-90985682447533284/44253206255743824031,
3311159023183359/6321886607963403433,
14020904561954201/6321886607963403433,
-30876860056242495/44253206255743824031,
1283188347895953/12643773215926806866,
5054950793447334/6321886607963403433,
-3714494491388754/6321886607963403433,
10363473683003709/44253206255743824031,
-44214967672541940/44253206255743824031,
3389336942373831/6321886607963403433,
-2097740039858619/12643773215926806866)
<b,y>:
-2097740039858619/12643773215926806866*2^(2/3)*3^(2/3)*5^(2/3) +
3311159023183359/6321886607963403433*2^(2/3)*3^(2/3)*5^(1/3) +
10363473683003709/44253206255743824031*2^(2/3)*3^(1/3)*5^(2/3) +
3389336942373831/6321886607963403433*2^(1/3)*3^(2/3)*5^(2/3) -
79218868344463653/88506412511487648062*2^(2/3)*3^(2/3) -
7004738263654053/12643773215926806866*2^(2/3)*3^(1/3)*5^(1/3) +
1283188347895953/12643773215926806866*2^(2/3)*5^(2/3) -
90985682447533284/44253206255743824031*2^(1/3)*3^(2/3)*5^(1/3) -
3714494491388754/6321886607963403433*2^(1/3)*3^(1/3)*5^(2/3) -
44214967672541940/44253206255743824031*3^(2/3)*5^(2/3) +
7609214068173315/12643773215926806866*2^(2/3)*3^(1/3) -
63303351032438961/88506412511487648062*2^(2/3)*5^(1/3) +
30694440328632370/6321886607963403433*2^(1/3)*3^(2/3) +
10778499819152451/6321886607963403433*2^(1/3)*3^(1/3)*5^(1/3) -
30876860056242495/44253206255743824031*2^(1/3)*5^(2/3) +
32291803150446358/6321886607963403433*3^(2/3)*5^(1/3) +
5054950793447334/6321886607963403433*3^(1/3)*5^(2/3) +
14159343943246004/6321886607963403433*2^(2/3) -
119414810568569283/44253206255743824031*2^(1/3)*3^(1/3) +
27999438968454556/6321886607963403433*2^(1/3)*5^(1/3) -
135774633227300964/44253206255743824031*3^(1/3)*5^(1/3) -
113737199574721883/6321886607963403433*2^(1/3) -
118070567437573862/6321886607963403433*3^(2/3) +
45903962039615943/6321886607963403433*3^(1/3) +
14020904561954201/6321886607963403433*5^(2/3) -
114594372661917362/6321886607963403433*5^(1/3) +
6085751511148292591/44253206255743824031

Basis Vector (b):
(1, 2^(1/3), 2^(2/3), 3^(1/3), 2^(1/3)*3^(1/3), 2^(2/3)*3^(1/3), 3^(2/3), 2^(1/3)*3^(2/3), 2^(2/3)*3^(2/3), 5^(1/3), 2^(1/3)*5^(1/3), 2^(2/3)*5^(1/3), 3^(1/3)*5^(1/3), 2^(1/3)*3^(1/3)*5^(1/3), 2^(2/3)*3^(1/3)*5^(1/3), 3^(2/3)*5^(1/3), 2^(1/3)*3^(2/3)*5^(1/3), 2^(2/3)*3^(2/3)*5^(1/3), 5^(2/3), 2^(1/3)*5^(2/3), 2^(2/3)*5^(2/3), 3^(1/3)*5^(2/3), 2^(1/3)*3^(1/3)*5^(2/3), 2^(2/3)*3^(1/3)*5^(2/3), 3^(2/3)*5^(2/3), 2^(1/3)*3^(2/3)*5^(2/3), 2^(2/3)*3^(2/3)*5^(2/3))
Expression vector (x):
(7, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
<x,b>:
2^(1/3) + 3^(2/3) + 5^(1/3) + 7
y:
(6085751511148292591/44253206255743824031, -113737199574721883/6321886607963403433, 14159343943246004/6321886607963403433, 45903962039615943/6321886607963403433, -119414810568569283/44253206255743824031, 7609214068173315/12643773215926806866, -118070567437573862/6321886607963403433, 30694440328632370/6321886607963403433, -79218868344463653/88506412511487648062, -114594372661917362/6321886607963403433, 27999438968454556/6321886607963403433, -63303351032438961/88506412511487648062, -135774633227300964/44253206255743824031, 10778499819152451/6321886607963403433, -7004738263654053/12643773215926806866, 32291803150446358/6321886607963403433, -90985682447533284/44253206255743824031, 3311159023183359/6321886607963403433, 14020904561954201/6321886607963403433, -30876860056242495/44253206255743824031, 1283188347895953/12643773215926806866, 5054950793447334/6321886607963403433, -3714494491388754/6321886607963403433, 10363473683003709/44253206255743824031, -44214967672541940/44253206255743824031, 3389336942373831/6321886607963403433, -2097740039858619/12643773215926806866)
<b,y>:
-2097740039858619/12643773215926806866*2^(2/3)*3^(2/3)*5^(2/3) + 3311159023183359/6321886607963403433*2^(2/3)*3^(2/3)*5^(1/3) + 10363473683003709/44253206255743824031*2^(2/3)*3^(1/3)*5^(2/3) + 3389336942373831/6321886607963403433*2^(1/3)*3^(2/3)*5^(2/3) - 79218868344463653/88506412511487648062*2^(2/3)*3^(2/3) - 7004738263654053/12643773215926806866*2^(2/3)*3^(1/3)*5^(1/3) + 1283188347895953/12643773215926806866*2^(2/3)*5^(2/3) - 90985682447533284/44253206255743824031*2^(1/3)*3^(2/3)*5^(1/3) - 3714494491388754/6321886607963403433*2^(1/3)*3^(1/3)*5^(2/3) - 44214967672541940/44253206255743824031*3^(2/3)*5^(2/3) + 7609214068173315/12643773215926806866*2^(2/3)*3^(1/3) - 63303351032438961/88506412511487648062*2^(2/3)*5^(1/3) + 30694440328632370/6321886607963403433*2^(1/3)*3^(2/3) + 10778499819152451/6321886607963403433*2^(1/3)*3^(1/3)*5^(1/3) - 30876860056242495/44253206255743824031*2^(1/3)*5^(2/3) + 32291803150446358/6321886607963403433*3^(2/3)*5^(1/3) + 5054950793447334/6321886607963403433*3^(1/3)*5^(2/3) + 14159343943246004/6321886607963403433*2^(2/3) - 119414810568569283/44253206255743824031*2^(1/3)*3^(1/3) + 27999438968454556/6321886607963403433*2^(1/3)*5^(1/3) - 135774633227300964/44253206255743824031*3^(1/3)*5^(1/3) - 113737199574721883/6321886607963403433*2^(1/3) - 118070567437573862/6321886607963403433*3^(2/3) + 45903962039615943/6321886607963403433*3^(1/3) + 14020904561954201/6321886607963403433*5^(2/3) - 114594372661917362/6321886607963403433*5^(1/3) + 6085751511148292591/44253206255743824031