"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/Ser.c" (14 Jan 2021, 7314 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Ser.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

    1 /* Copyright (C) 2000, 2012  The PARI group.
    2 
    3 This file is part of the PARI/GP package.
    4 
    5 PARI/GP is free software; you can redistribute it and/or modify it under the
    6 terms of the GNU General Public License as published by the Free Software
    7 Foundation; either version 2 of the License, or (at your option) any later
    8 version. It is distributed in the hope that it will be useful, but WITHOUT
    9 ANY WARRANTY WHATSOEVER.
   10 
   11 Check the License for details. You should have received a copy of it, along
   12 with the package; see the file 'COPYING'. If not, write to the Free Software
   13 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
   14 
   15 #include "pari.h"
   16 #include "paripriv.h"
   17 /*******************************************************************/
   18 /*                                                                 */
   19 /*                     Conversion --> t_SER                        */
   20 /*                                                                 */
   21 /*******************************************************************/
   22 static GEN
   23 RgX_to_ser_i(GEN x, long l, long v, int copy)
   24 {
   25   long i = 2, lx = lg(x), vx = varn(x);
   26   GEN y;
   27   if (lx == 2) return zeroser(vx, minss(l - 2, v));
   28   if (l <= 2)
   29   {
   30     if (l == 2 && v != LONG_MAX) return zeroser(vx, v);
   31     pari_err_BUG("RgX_to_ser (l < 2)");
   32   }
   33   y = cgetg(l,t_SER);
   34   if (v == LONG_MAX) { v = 1; lx = 3; } /* e.g. Mod(0,3) * x^0 */
   35   else if (v)
   36   {
   37     long w = v;
   38     while (isrationalzero(gel(x,2))) { x++; w--; }
   39     lx -= v;
   40     if (w)
   41     { /* keep type information, e.g. Mod(0,3) + x */
   42       GEN z = gel(x,2); /* = 0 */
   43       x += w; gel(y,2) = gadd(gel(x,2), z); i++;
   44     }
   45   }
   46   y[1] = evalvarn(vx) | evalvalp(v); /* must come here because of LONG_MAX */
   47   if (lx > l) lx = l;
   48   /* 3 <= lx <= l */
   49   if (copy)
   50     for (; i <lx; i++) gel(y,i) = gcopy(gel(x,i));
   51   else
   52     for (; i <lx; i++) gel(y,i) = gel(x,i);
   53   for (     ; i < l; i++) gel(y,i) = gen_0;
   54   return normalize(y);
   55 }
   56 /* enlarge/truncate t_POL x to a t_SER with lg l */
   57 GEN
   58 RgX_to_ser(GEN x, long l) { return RgX_to_ser_i(x, l, RgX_val(x), 0); }
   59 GEN
   60 RgX_to_ser_inexact(GEN x, long l)
   61 {
   62   long i, lx = lg(x);
   63   int first = 1;
   64   for (i = 2; i < lx && gequal0(gel(x,i)); i++) /* ~ RgX_valrem + normalize */
   65     if (first && !isexactzero(gel(x,i)))
   66     {
   67       pari_warn(warner,"normalizing a series with 0 leading term");
   68       first = 0;
   69     }
   70   return RgX_to_ser_i(x, l, i - 2, 0);
   71 }
   72 GEN
   73 rfrac_to_ser(GEN x, long l)
   74 {
   75   GEN d = gel(x,2);
   76   if (l == 2)
   77   {
   78     long v = varn(d);
   79     return zeroser(varn(d), gvaluation(x, pol_x(v)));
   80   }
   81   return gdiv(gel(x,1), RgX_to_ser(d, l));
   82 }
   83 
   84 static GEN
   85 RgV_to_ser_i(GEN x, long v, long l, int copy)
   86 {
   87   long j, lx = minss(lg(x), l-1);
   88   GEN y;
   89   if (lx == 1) return zeroser(v, l-2);
   90   y = cgetg(l, t_SER); y[1] = evalvarn(v)|evalvalp(0);
   91   x--;
   92   if (copy)
   93     for (j = 2; j <= lx; j++) gel(y,j) = gcopy(gel(x,j));
   94   else
   95     for (j = 2; j <= lx; j++) gel(y,j) = gel(x,j);
   96   for (     ; j < l;   j++) gel(y,j) = gen_0;
   97   return normalize(y);
   98 }
   99 GEN
  100 RgV_to_ser(GEN x, long v, long l) { return RgV_to_ser_i(x, v, l, 0); }
  101 
  102 /* x a t_SER, prec >= 0 */
  103 GEN
  104 sertoser(GEN x, long prec)
  105 {
  106   long i, lx = lg(x), l;
  107   GEN y;
  108   if (lx == 2) return zeroser(varn(x), prec);
  109   l = prec+2; lx = minss(lx, l);
  110   y = cgetg(l,t_SER); y[1] = x[1];
  111   for (i = 2; i < lx; i++) gel(y,i) = gel(x,i);
  112   for (     ; i < l;  i++) gel(y,i) = gen_0;
  113   return y;
  114 }
  115 
  116 /* R(1/x) = x^v * n/d, val(n) = val(d) = 0 */
  117 long
  118 rfracrecip(GEN *pn, GEN *pd)
  119 {
  120   long v = degpol(*pd);
  121   if (typ(*pn) == t_POL && varn(*pn) == varn(*pd))
  122   {
  123     v -= degpol(*pn);
  124     (void)RgX_valrem(*pn, pn); *pn = RgX_recip(*pn);
  125   }
  126   (void)RgX_valrem(*pd, pd); *pd = RgX_recip(*pd);
  127   return v;
  128 }
  129 
  130 /* R(1/x) + O(x^N) */
  131 GEN
  132 rfracrecip_to_ser_absolute(GEN R, long N)
  133 {
  134   GEN n = gel(R,1), d = gel(R,2);
  135   long v = rfracrecip(&n, &d); /* R(1/x) = x^v * n/d, val(n) = val(d) = 0 */
  136   if (N <= v) return zeroser(varn(d), N);
  137   R = gdiv(n, RgX_to_ser(d, N-v+2));
  138   setvalp(R, v); return R;
  139 }
  140 
  141 /* assume prec >= 0 */
  142 GEN
  143 scalarser(GEN x, long v, long prec)
  144 {
  145   long i, l;
  146   GEN y;
  147 
  148   if (gequal0(x))
  149   {
  150     if (isrationalzero(x)) return zeroser(v, prec);
  151     if (!isexactzero(x)) prec--;
  152     y = cgetg(3, t_SER);
  153     y[1] = evalsigne(0) | _evalvalp(prec) | evalvarn(v);
  154     gel(y,2) = gcopy(x); return y;
  155   }
  156   l = prec + 2; y = cgetg(l, t_SER);
  157   y[1] = evalsigne(1) | _evalvalp(0) | evalvarn(v);
  158   gel(y,2) = gcopy(x); for (i=3; i<l; i++) gel(y,i) = gen_0;
  159   return y;
  160 }
  161 
  162 /* assume x a t_[SER|POL], apply gtoser to all coeffs */
  163 static GEN
  164 coefstoser(GEN x, long v, long prec)
  165 {
  166   long i, lx;
  167   GEN y = cgetg_copy(x, &lx); y[1] = x[1];
  168   for (i=2; i<lx; i++) gel(y,i) = gtoser(gel(x,i), v, prec);
  169   return y;
  170 }
  171 
  172 static void
  173 err_ser_priority(GEN x, long v) { pari_err_PRIORITY("Ser", x, "<", v); }
  174 /* x a t_POL */
  175 static GEN
  176 poltoser(GEN x, long v, long prec)
  177 {
  178   long s = varncmp(varn(x), v);
  179   if (s < 0) err_ser_priority(x,v);
  180   if (s > 0) return scalarser(x, v, prec);
  181   return RgX_to_ser_i(x, prec+2, RgX_val(x), 1);
  182 }
  183 /* x a t_RFRAC */
  184 static GEN
  185 rfractoser(GEN x, long v, long prec)
  186 {
  187   long s = varncmp(varn(gel(x,2)), v);
  188   pari_sp av;
  189   if (s < 0) err_ser_priority(x,v);
  190   if (s > 0) return scalarser(x, v, prec);
  191   av = avma; return gerepileupto(av, rfrac_to_ser(x, prec+2));
  192 }
  193 GEN
  194 toser_i(GEN x)
  195 {
  196   switch(typ(x))
  197   {
  198     case t_SER: return x;
  199     case t_POL: return RgX_to_ser(x, precdl+2);
  200     case t_RFRAC: return rfrac_to_ser(x, precdl+2);
  201   }
  202   return NULL;
  203 }
  204 
  205 /* conversion: prec ignored if t_VEC or t_SER in variable v */
  206 GEN
  207 gtoser(GEN x, long v, long prec)
  208 {
  209   long tx = typ(x);
  210 
  211   if (v < 0) v = 0;
  212   if (prec < 0) pari_err_DOMAIN("Ser", "precision", "<", gen_0, stoi(prec));
  213   if (tx == t_SER)
  214   {
  215     long s = varncmp(varn(x), v);
  216     if      (s < 0) return coefstoser(x, v, prec);
  217     else if (s > 0) return scalarser(x, v, prec);
  218     return gcopy(x);
  219   }
  220   if (is_scalar_t(tx)) return scalarser(x,v,prec);
  221   switch(tx)
  222   {
  223     case t_POL: return poltoser(x, v, prec);
  224     case t_RFRAC: return rfractoser(x, v, prec);
  225     case t_QFR:
  226     case t_QFI: return RgV_to_ser_i(x, v, 4+1, 1);
  227     case t_VECSMALL: x = zv_to_ZV(x);/*fall through*/
  228     case t_VEC: case t_COL:
  229       if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("Ser", x, "<=", v);
  230       return RgV_to_ser_i(x, v, lg(x)+1, 1);
  231   }
  232   pari_err_TYPE("Ser",x);
  233   return NULL; /* LCOV_EXCL_LINE */
  234 }
  235 /* impose prec */
  236 GEN
  237 gtoser_prec(GEN x, long v, long prec)
  238 {
  239   pari_sp av = avma;
  240   if (v < 0) v = 0;
  241   if (prec < 0) pari_err_DOMAIN("Ser", "precision", "<", gen_0, stoi(prec));
  242   switch(typ(x))
  243   {
  244     case t_SER: if (varn(x) != v) break;
  245                 return gerepilecopy(av, sertoser(x, prec));
  246     case t_QFR:
  247     case t_QFI:
  248       x = RgV_to_ser_i(mkvec3(gel(x,1),gel(x,2),gel(x,3)), v, prec+2, 1);
  249       return gerepileupto(av, x);
  250     case t_VECSMALL: x = zv_to_ZV(x);/*fall through*/
  251     case t_VEC: case t_COL:
  252       if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("Ser", x, "<=", v);
  253       return RgV_to_ser_i(x, v, prec+2, 1);
  254   }
  255   return gtoser(x, v, prec);
  256 }
  257 GEN
  258 Ser0(GEN x, long v, GEN d, long prec)
  259 {
  260   if (!d) return gtoser(x, v, prec);
  261   if (typ(d) != t_INT)
  262   {
  263     d = gceil(d);
  264     if (typ(d) != t_INT) pari_err_TYPE("Ser [precision]",d);
  265   }
  266   return gtoser_prec(x, v, itos(d));
  267 }