share/pslq/pslq.mac

   1 /***********************************
   2  **
   3  ** Integer relations package (pslq.mac)
   4  **
   5  ** Copyright: 2006, Andrej Vodopivec <andrej.vodopivec@gmail.com>
   6  ** Licence:   GPL
   7  **
   8  **   pslq.mac implements guess_exact_value function, which is used
   9  **   to identify floating point numbers as functions of known constatns.
  10  **
  11  ** For detailed description see
  12  **
  13  **  P.Borwein, K.G.Hare, A.Meichsner: Reverse symbolic computations, the
  14  **     identify function.
  15  **
  16  ***************************
  17   Examples:
  18
  19   (%i2) float(%e^2+%e-1);
  20   (%o2) 9.107337927389695
  21   (%i3) guess_exact_value(%);
  22   (%o3) %e^2+%e-1
  23   (%i4) float(sin(%pi/30));
  24   (%o4) 0.10452846326765
  25   (%i5) guess_exact_value(%);
  26   (%o5) (sqrt(sqrt(5)-1)*sqrt(15))/(4*sqrt(2)*5^(1/4))-sqrt(5)/8-1/8
  27   (%i6) float(%pi*%e/2);
  28   (%o6) 4.269867111336783
  29   (%i7) guess_exact_value(%);
  30   (%o7) (%e*%pi)/2
  31
  32  ***********************/
  33
  34 define_variable(identify_sum_constants,
  35   [%pi, %e, sqrt(2), sqrt(3), sqrt(5), log(2), log(3)], list)$
  36 define_variable(identify_product_constants,
  37   [%pi, %e, sqrt(2), sqrt(3), sqrt(5), log(2), log(3)], list)$
  38 define_variable(identify_powersum_constants,
  39   [%pi, %e, log(2), log(3)], list)$
  40 define_variable(identify_extend,
  41   [], list)$
  42 define_variable(identify_powersum_max_degree,
  43   3, fixnum)$
  44 define_variable(identify_powersum_min_degree,
  45   2, fixnum)$
  46 define_variable(identify_algebraic_degree,
  47   4, fixnum)$
  48
  49 if get(pslq_integer_relation, version)=false
  50   then if load_pathname # false
  51          then load (sconcat (pathname_directory (load_pathname), "pslq_integer_relation.lisp"))
  52        else load ("pslq_integer_relation.lisp");
  53
  54 identify_integer_relation(x, [threshold]) := block(
  55   [relation : pslq_integer_relation(x)],
  56   if length(threshold)#1 then threshold:100 else threshold:threshold[1],
  57   if relation=false then return(false),
  58   if apply("+", (map(abs, relation)))>threshold then return(false),
  59   relation
  60 )$
  61
  62 guess_exact_value(x) :=
  63   if mapatom(x) then (
  64     if floatnump(x) or bfloatp(x) then block(
  65       [res],
  66       if equal(x, 0) or equal(x, -1) or equal(x, 1)
  67         then return(floor(x)),
  68       if errcatch (
  69         res : identify_as_algebraic(x)
  70       ) = [] then res : false,
  71       if res#false then return(res),
  72       if errcatch (
  73         res : identify_as_sum(x)
  74       ) = [] then res : false,
  75       if res#false then return(res),
  76       if errcatch (
  77         res : identify_as_product(x)
  78       ) = [] then res : false,
  79       if res#false then return(res),
  80       if errcatch (
  81         res : identify_as_powersum(x)
  82       ) = [] then res : false,
  83       if res#false then return(res),
  84       x
  85     )
  86     else x
  87   )
  88   else
  89     map(guess_exact_value, x)$
  90
  91 identify_as_sum(x) := block(
  92   [floats, relation, res : false, i,
  93   identify_sum_constants : append(identify_sum_constants, identify_extend)],
  94   if floatnump(x) then floats : append([x], float(identify_sum_constants))
  95   else (
  96     x : bfloat(x),
  97     if not(bfloatp(x)) then error("Argument to `identify_as_sum' is not a float!"),
  98     floats : append([x], bfloat(identify_sum_constants))
  99   ),
 100   relation : identify_integer_relation(floats),
 101
 102   if relation#false then (
 103     res : 0,
 104     for i:1 thru length(identify_sum_constants) do (
 105       res : res + relation[i+1]*identify_sum_constants[i]
 106     ),
 107     res : -res/relation[1]
 108   ),
 109
 110   res
 111 )$
 112
 113 identify_as_algebraic(x) := block(
 114   [floats, relation, res : false, i, pol, %x%, dif, r, rr, ratprint:false, sol],
 115   if floatnump(x) then floats : makelist(x^i, i, 0, identify_algebraic_degree)
 116   else (
 117     x : bfloat(x),
 118     if not(bfloatp(x)) then error("Argument to `identify_as_algebraic' is not a float!"),
 119     floats : makelist(x^i, i, 0, identify_algebraic_degree)
 120   ),
 121   relation : identify_integer_relation(floats, 100*identify_algebraic_degree),
 122
 123   if relation#false then (
 124     pol : makelist(%x%^i, i, 0, identify_algebraic_degree) . relation,
 125     sol : solve(pol, %x%),
 126     for r in sol do (
 127       r : rhs(r),
 128       if numberp(rectform(float(r))) then (
 129         if floatnump(x) then rr : float(r) else rr : bfloat(r),
 130         if not(floatnump(r)) and cabs(rr-x)<10^-14 then res : r
 131       )
 132     )
 133   ),
 134
 135   res
 136 )$
 137
 138 identify_as_product(x) := block(
 139   [floats, relation, res : false,
 140   identify_product_constants : append(identify_product_constants, identify_extend)],
 141   if floatnump(x) then floats : append([log(abs(x))], float(map(log, identify_product_constants)))
 142   else (
 143     x : bfloat(x),
 144     if not(bfloatp(x)) then error("Argument to `identify_as_sum' is not a float!"),
 145     floats : append([log(abs(x))], bfloat(map(log, identify_product_constants)))
 146   ),
 147   relation : identify_integer_relation(floats),
 148
 149   if relation#false then (
 150     res : 1,
 151     for i:1 thru length(identify_product_constants) do (
 152       res : res * identify_product_constants[i]^relation[i+1]
 153       ),
 154     res : res^(-1/relation[1])
 155   )
 156   else return(false),
 157
 158   if x<0 then
 159     -res
 160   else
 161     res
 162 )$
 163
 164 identify_as_powersum(x) := block(
 165   [res : false, floats, relation,
 166   identify_powersum_constants : append(identify_powersum_constants, identify_extend),
 167   min_deg : identify_powersum_min_degree,
 168   max_deg : identify_powersum_max_degree],
 169   if not(floatnump(x)) then (
 170     x : bfloat(x),
 171     if not(bfloatp(x)) then error("Argument to `identify_as_powersum' is not a float!")
 172   ),
 173
 174   res : for c in identify_powersum_constants do (
 175     if floatnump(x) then
 176       floats : makelist(float(c^i), i, -min_deg, max_deg)
 177     else
 178       floats : makelist(bfloat(c^i), i, -min_deg, max_deg),
 179     floats : append([x], floats),
 180     relation : identify_integer_relation(floats),
 181     if relation#false and relation[1]#0 then (
 182       res : makelist(-relation[i+1]/relation[1], i,
 183         1, min_deg + max_deg + 1),
 184       res : res . makelist(c^i, i, -min_deg, max_deg),
 185       return(res)
 186     )
 187   ),
 188
 189   if res=done then return(false),
 190   res
 191 )$