bridge/class-master/depend/pari/src/language/sumiter.c
/* Find zeros of a function in the real interval [a,b] by interval splitting */
Find zeros of a function in the real interval [a,b] by interval splitting
Find zeros of a function in the real interval [ a b ] by interval splitting
GEN
solvestep(void *E, GEN (*f)(void *,GEN), GEN a, GEN b, GEN step, long flag, long prec)
{
  const long ITMAX = 10;
  pari_sp av = avma;
  GEN fa, ainit, binit;
  long sainit, it, bit = prec2nbits(prec) / 2, ct = 0, s = gcmp(a,b);
  if (!s) return gequal0(f(E, a)) ? gcopy(mkvec(a)): cgetg(1,t_VEC);
  if (s > 0) swap(a, b);
  if (flag&4)
  {
    if (gcmpgs(step,1)<=0) pari_err_DOMAIN("solvestep","step","<=",gen_1,step);
    if (gsigne(a) <= 0) pari_err_DOMAIN("solvestep","a","<=",gen_0,a);
  }
  else if (gsigne(step) <= 0)
    pari_err_DOMAIN("solvestep","step","<=",gen_0,step);
  ainit = a = gtofp(a, prec); fa = f(E, a);
  binit = b = gtofp(b, prec); step = gtofp(step, prec);
  sainit = gsigne(fa);
  if (gexpo(fa) < -bit) sainit = 0;
  for (it = 0; it < ITMAX; it++)
  {
    pari_sp av2 = avma;
    GEN v = cgetg(1, t_VEC);
    long sa;
    a = ainit;
    b = binit;
    sa = sainit;
    while (gcmp(a,b) < 0)
    {
      GEN fc, c = (flag&4)? gmul(a, step): gadd(a, step);
      long sc;
      if (gcmp(c,b) > 0) c = b;
      fc = f(E, c);
      sc = gsigne(fc);
      if (gexpo(fc) < -bit) sc = 0;
      if (!sc || sa*sc < 0)
      {
        long e;
        GEN z;
        z = sc? zbrent(E, f, a, c, prec): c;
        (void)grndtoi(z, &e);
        if (e  <= -bit) ct = 1;
        if ((flag&1) && ((!(flag&8)) || ct)) return gerepileupto(av, z);
        v = gconcat(v, z);
      }
      a = c; fa = fc; sa = sc;
    }
    if ((!(flag&2) || lg(v) > 1) && (!(flag&8) || ct))
      return gerepilecopy(av, v);
    step = (flag&4)? sqrtr(sqrtr(step)): gmul2n(step, -2);
    gerepileall(av2, 2, &fa, &step);
  }
  if (it == ITMAX) pari_err_IMPL("solvestep recovery [too many iterations]");
  return NULL;
}
GEN solvestep ( void * E , GEN ( * f ) ( void * , GEN ) , GEN a , GEN b , GEN step , long flag , long prec ) { const long ITMAX = 10 ; pari_sp av = avma ; GEN fa , ainit , binit ; long sainit , it , bit = prec2nbits ( prec ) / 2 , ct = 0 , s = gcmp ( a , b ) ; if ( ! s ) return gequal0 ( f ( E , a ) ) ? gcopy ( mkvec ( a ) ) : cgetg ( 1 , t_VEC ) ; if ( s > 0 ) swap ( a , b ) ; if ( flag & 4 ) { if ( gcmpgs ( step , 1 ) <= 0 ) pari_err_DOMAIN ( " " , " " , " " , gen_1 , step ) ; if ( gsigne ( a ) <= 0 ) pari_err_DOMAIN ( " " , " " , " " , gen_0 , a ) ; } else if ( gsigne ( step ) <= 0 ) pari_err_DOMAIN ( " " , " " , " " , gen_0 , step ) ; ainit = a = gtofp ( a , prec ) ; fa = f ( E , a ) ; binit = b = gtofp ( b , prec ) ; step = gtofp ( step , prec ) ; sainit = gsigne ( fa ) ; if ( gexpo ( fa ) < - bit ) sainit = 0 ; for ( it = 0 ; it < ITMAX ; it ++ ) { pari_sp av2 = avma ; GEN v = cgetg ( 1 , t_VEC ) ; long sa ; a = ainit ; b = binit ; sa = sainit ; while ( gcmp ( a , b ) < 0 ) { GEN fc , c = ( flag & 4 ) ? gmul ( a , step ) : gadd ( a , step ) ; long sc ; if ( gcmp ( c , b ) > 0 ) c = b ; fc = f ( E , c ) ; sc = gsigne ( fc ) ; if ( gexpo ( fc ) < - bit ) sc = 0 ; if ( ! sc || sa * sc < 0 ) { long e ; GEN z ; z = sc ? zbrent ( E , f , a , c , prec ) : c ; ( void ) grndtoi ( z , & e ) ; if ( e <= - bit ) ct = 1 ; if ( ( flag & 1 ) && ( ( ! ( flag & 8 ) ) || ct ) ) return gerepileupto ( av , z ) ; v = gconcat ( v , z ) ; } a = c ; fa = fc ; sa = sc ; } if ( ( ! ( flag & 2 ) || lg ( v ) > 1 ) && ( ! ( flag & 8 ) || ct ) ) return gerepilecopy ( av , v ) ; step = ( flag & 4 ) ? sqrtr ( sqrtr ( step ) ) : gmul2n ( step , -2 ) ; gerepileall ( av2 , 2 , & fa , & step ) ; } if ( it == ITMAX ) pari_err_IMPL ( " " ) ; return NULL ; }
Find zeros of a function in the real interval [a,b] by interval splitting
Find zeros of a function in the real interval [ a b ] by interval splitting