147 lines
4.9 KiB
Plaintext
147 lines
4.9 KiB
Plaintext
/*
|
|
Copyright (c) 2015 Thomas Baruchel
|
|
|
|
Permission is hereby granted, free of charge, to any person obtaining a copy of
|
|
this software and associated documentation files (the "Software"), to deal in
|
|
the Software without restriction, including without limitation the rights to
|
|
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
|
|
of the Software, and to permit persons to whom the Software is furnished to do
|
|
so, subject to the following conditions:
|
|
|
|
The above copyright notice and this permission notice shall be included in all
|
|
copies or substantial portions of the Software.
|
|
|
|
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
SOFTWARE.
|
|
*/
|
|
|
|
/*
|
|
|
|
An integer relation algorithm for finding integer relations between constants
|
|
(bfloat type being the intended type to be used by the function).
|
|
|
|
See: https://en.wikipedia.org/wiki/Integer_relation_algorithm
|
|
|
|
Examples
|
|
========
|
|
|
|
(%i4) fpprec: 64$
|
|
(%i5) pslq([bfloat(zeta(2)), %pi^2], 1b-32);
|
|
(%o5) [- 6, 1]
|
|
|
|
(%i78) fpprec: 64$ phi: bfloat((1+sqrt(5))/2)$
|
|
(%i79) algdep(phi,2,1b-32);
|
|
(%o79) x^2-x-1
|
|
|
|
(%i65) fpprec: 64$ zerobern: false$
|
|
(%i67) b:makelist(bfloat(bern(k)), k, 10);
|
|
(%o67) [- 5.0b-1,
|
|
1.666666666666666666666666666666666666666666666666666666666666667b-1,
|
|
- 3.333333333333333333333333333333333333333333333333333333333333333b-2,
|
|
7.575757575757575757575757575757575757575757575757575757575757576b-2,
|
|
- 2.531135531135531135531135531135531135531135531135531135531135531b-1,
|
|
1.166666666666666666666666666666666666666666666666666666666666667b0,
|
|
- 7.092156862745098039215686274509803921568627450980392156862745098b0,
|
|
5.497117794486215538847117794486215538847117794486215538847117795b1,
|
|
- 5.291242424242424242424242424242424242424242424242424242424242424b2,
|
|
6.192123188405797101449275362318840579710144927536231884057971015b3]
|
|
(%i68) makelist( algdep( b[k], 1, 1b-32), k, 10);
|
|
(%o68) [2 x + 1, 1 - 6 x, 30 x + 1, 66 x - 5, - 2730 x - 691, 7 - 6 x,
|
|
510 x + 3617, 43867 - 798 x, 330 x + 174611, 138 x - 854513]
|
|
*/
|
|
|
|
define_variable(pslqMaxItr, 256, integer)$
|
|
|
|
/* TODO:
|
|
known bugs:
|
|
(1) fpprec: 64$
|
|
pslq([bfloat(1/90),1],1b-32); Works
|
|
pslq([1,bfloat(1/90)],1b-32); Fails
|
|
|
|
|
|
*/
|
|
pslq(x, prec) := block([
|
|
n:length(x),
|
|
teps:prec*1.6b1,
|
|
gam: 1.2b0,
|
|
t, y, s, h, a, b, m,
|
|
done: false
|
|
],
|
|
/* Compute the vector s */
|
|
s:makelist(0.0b0, n),
|
|
s[n]: bfloat(abs( x[n] )),
|
|
t: x[n]*x[n],
|
|
for i:(n-1) step -1 thru 1 do (
|
|
t: bfloat(t + x[i]*x[i]),
|
|
s[i]: sqrt(t) ),
|
|
/* Normalize the vector x, s */
|
|
t: s[1],
|
|
y: makelist( bfloat(x[i]/t), i, n),
|
|
s: map(lambda([z], z/t), s),
|
|
/* Construct matrix H, a, b */
|
|
h: zeromatrix(n, n-1),
|
|
a: ident(n),
|
|
b: ident(n),
|
|
for i:1 thru n do
|
|
for j:1 thru min(i, n-1) do
|
|
h[i,j]: if i=j then s[j+1]/s[j] else -y[i]*y[j]/(s[j]*s[j+1]),
|
|
/* Reduce matrix H */
|
|
for i:2 thru n do
|
|
for j:(i-1) step -1 thru 1 do (
|
|
t: round( h[i,j]/h[j,j] ),
|
|
y[j]: y[j]+t*y[i],
|
|
for k:1 thru j do h[i,k]: h[i,k] - t *h[j,k],
|
|
for k:1 thru n do (
|
|
a[i,k]: a[i,k] - t * a[j,k],
|
|
b[k,j]: b[k,j] + t * b[k,i] ) ),
|
|
|
|
for itr:1 while (itr <= pslqMaxItr) and (not done)
|
|
do block([ mval:-1, g:gam ],
|
|
for i:1 thru n-1 do (
|
|
t: abs(g * h[i,i]),
|
|
if t > mval then ( mval:t, m:i ),
|
|
g: g*gam),
|
|
t: y[m], y[m]: y[m+1], y[m+1]: t,
|
|
a:rowswap(a, m, m+1),
|
|
h:rowswap(h, m, m+1),
|
|
b:columnswap(b,m, m+1),
|
|
if m < n-1 then block([t0: sqrt( h[m,m]^2 + h[m,m+1]^2), t1, t2, t3, t4],
|
|
t1: h[m,m]/ t0,
|
|
t2: h[m,m+1]/ t0,
|
|
for i:m thru n do (
|
|
t3: h[i,m],
|
|
t4: h[i,m+1],
|
|
h[i,m]: t1*t3 + t2*t4,
|
|
h[i,m+1]: t1*t4 - t2*t3)
|
|
),
|
|
for i:m+1 thru n do
|
|
for j:min(i-1,m+1) step -1 thru 1 do (
|
|
t: round( h[i,j]/h[j,j] ),
|
|
y[j]: y[j] + t * y[i],
|
|
for k:1 thru j do h[i,k]: h[i,k] - t*h[j,k],
|
|
for k:1 thru n do (
|
|
a[i,k]: a[i,k] - t * a[j,k],
|
|
b[k,j]: b[k,j] + t * b[k,i] )),
|
|
|
|
mval: 0,
|
|
for j: 1 thru n-1 do (
|
|
t: abs( h[j,j] ),
|
|
if t > mval then mval: t),
|
|
for i:1 thru n do (
|
|
t: abs( y[i] ),
|
|
if t < teps then (done: true, m: i))
|
|
),
|
|
if done then makelist(b[i,m],i,n) else makelist(0, n)
|
|
)$
|
|
|
|
algdep(x, n, prec) := block([k, v],
|
|
v: makelist(1, n+1), for k:n thru 1 step -1 do v[k]:v[k+1]*x,
|
|
v: pslq(v, prec),
|
|
sum( v[i] * ('x)^(n+1-i), i, 1, n+1)
|
|
)$
|