Ticket #13899: 13899_TAB_sagelib.patch

File 13899_TAB_sagelib.patch, 148.1 KB (added by jdemeyer, 7 years ago)
  • c_lib/include/ZZ_pylong.h

    # HG changeset patch
    # User Jeroen Demeyer <jdemeyer@cage.ugent.be>
    # Date 1357126969 -3600
    # Node ID b57bdc0626bf007de4720079553c9b3ae453d8c2
    # Parent  da56551f7e8fe3245ff30480f65be6eadc327564
    Fix indentation and trailing spaces
    
    diff --git a/c_lib/include/ZZ_pylong.h b/c_lib/include/ZZ_pylong.h
    a b  
    1515#                  http://www.gnu.org/licenses/
    1616#*****************************************************************************/
    1717
    18 /*  Author:  Joel B. Mohler <joel@kiwistrawberry.us>
    19                          2007-06-17                                                      */
     18/*  Author:  Joel B. Mohler <joel@kiwistrawberry.us> (2007-06-17) */
    2019
    2120#ifndef ZZ_PYLONG_H
    2221#define ZZ_PYLONG_H
  • c_lib/include/ginac_wrap.h

    diff --git a/c_lib/include/ginac_wrap.h b/c_lib/include/ginac_wrap.h
    a b  
    6363
    6464bool relational_to_bool(const ex& e) {
    6565    if (ex_to<relational>(e))
    66         return 1;
     66        return 1;
    6767    else
    68         return 0;
     68        return 0;
    6969}
    7070
    7171bool g_is_a_terminating_series(const ex& e) {
    7272    if (is_a<pseries>(e)) {
    73         return (ex_to<pseries>(e)).is_terminating();
     73        return (ex_to<pseries>(e)).is_terminating();
    7474    }
    7575    return false;
    7676}
     
    8383relational::operators switch_operator(relational::operators o) {
    8484    switch(o) {
    8585    case relational::less:
    86         return relational::greater;
     86        return relational::greater;
    8787    case relational::less_or_equal:
    88         return relational::greater_or_equal;
     88        return relational::greater_or_equal;
    8989    case relational::greater:
    90         return relational::less;
     90        return relational::less;
    9191    case relational::greater_or_equal:
    92         return relational::less_or_equal;
     92        return relational::less_or_equal;
    9393    default:
    94         return o;
     94        return o;
    9595    }
    9696}
    9797
    9898bool is_negative(ex x) {
    9999    if (is_a<numeric>(x)) {
    100         return (ex_to<numeric>(x)).is_negative();
     100        return (ex_to<numeric>(x)).is_negative();
    101101    }
    102102    return false;
    103103}
  • c_lib/src/ZZ_pylong.cpp

    diff --git a/c_lib/src/ZZ_pylong.cpp b/c_lib/src/ZZ_pylong.cpp
    a b  
    77#                   http://www.gnu.org/licenses/
    88#*****************************************************************************/
    99
    10 /*  Author:  Joel B. Mohler <joel@kiwistrawberry.us>
    11                          2007-06-17                                                      */
     10/*  Author:  Joel B. Mohler <joel@kiwistrawberry.us> (2007-06-17) */
    1211
    1312#include "ZZ_pylong.h"
    1413#include "ntl_wrap.h"
  • c_lib/src/ntl_wrap.cpp

    diff --git a/c_lib/src/ntl_wrap.cpp b/c_lib/src/ntl_wrap.cpp
    a b  
    88#include <NTL/tools.h>
    99
    1010void del_charstar(char* a) {
    11   delete[] a;
     11    delete[] a;
    1212}
    1313
    1414
    1515void setup_NTL_error_callback(void (*function)(const char*, void*), void* context)
    1616{
    17    NTL::SetErrorCallbackFunction(function, context);
     17    NTL::SetErrorCallbackFunction(function, context);
    1818}
    1919
    2020
     
    2424   AUTHOR: David Harvey (2008-06-08) */
    2525int ZZ_to_int(const ZZ* x)
    2626{
    27   return to_int(*x);
     27    return to_int(*x);
    2828}
    2929
    3030/* Returns a *new* ZZ object.
    3131   AUTHOR: David Harvey (2008-06-08) */
    3232struct ZZ* int_to_ZZ(int value)
    3333{
    34   ZZ* output = new ZZ();
    35   conv(*output, value);
    36   return output;
     34    ZZ* output = new ZZ();
     35    conv(*output, value);
     36    return output;
    3737}
    3838
    3939/* Copies the ZZ into the mpz_t
     
    8383   AUTHOR: David Harvey (2008-06-08) */
    8484void ZZ_set_from_int(ZZ* x, int value)
    8585{
    86   conv(*x, value);
     86    conv(*x, value);
    8787}
    8888
    8989long ZZ_remove(struct ZZ &dest, const struct ZZ &src, const struct ZZ &f)
     
    153153   AUTHOR: David Harvey (2008-06-08) */
    154154int ZZ_p_to_int(const ZZ_p& x )
    155155{
    156   return ZZ_to_int(&rep(x));
     156    return ZZ_to_int(&rep(x));
    157157}
    158158
    159159/* Returns a *new* ZZ_p object.
    160160   AUTHOR: David Harvey (2008-06-08) */
    161161ZZ_p int_to_ZZ_p(int value)
    162162{
    163   ZZ_p r;
    164   r = value;
    165   return r;
     163    ZZ_p r;
     164    r = value;
     165    return r;
    166166}
    167167
    168168/* Sets given ZZ_p to value
    169169   AUTHOR: David Harvey (2008-06-08) */
    170170void ZZ_p_set_from_int(ZZ_p* x, int value)
    171171{
    172   conv(*x, value);
     172    conv(*x, value);
    173173}
    174174
    175175void ZZ_p_modulus(struct ZZ* mod, const struct ZZ_p* x)
    176176{
    177         (*mod) = x->modulus();
     177    (*mod) = x->modulus();
    178178}
    179179
    180180struct ZZ_p* ZZ_p_pow(const struct ZZ_p* x, long e)
    181181{
    182   ZZ_p *z = new ZZ_p();
    183   power(*z, *x, e);
    184   return z;
     182    ZZ_p *z = new ZZ_p();
     183    power(*z, *x, e);
     184    return z;
    185185}
    186186
    187187void ntl_ZZ_set_modulus(ZZ* x)
    188188{
    189   ZZ_p::init(*x);
     189    ZZ_p::init(*x);
    190190}
    191191
    192192ZZ_p* ZZ_p_inv(ZZ_p* x)
    193193{
    194   ZZ_p *z = new ZZ_p();
    195   inv(*z, *x);
    196   return z;
     194    ZZ_p *z = new ZZ_p();
     195    inv(*z, *x);
     196    return z;
    197197}
    198198
    199199ZZ_p* ZZ_p_random(void)
    200200{
    201   ZZ_p *z = new ZZ_p();
    202   random(*z);
    203   return z;
     201    ZZ_p *z = new ZZ_p();
     202    random(*z);
     203    return z;
    204204}
    205205
    206206struct ZZ_p* ZZ_p_neg(struct ZZ_p* x)
    207207{
    208   return new ZZ_p(-(*x));
     208    return new ZZ_p(-(*x));
    209209}
    210210
    211211
     
    216216
    217217char* ZZX_repr(struct ZZX* x)
    218218{
    219   ostringstream instore;
    220   instore << (*x);
    221   int n = strlen(instore.str().data());
    222   char* buf = new char[n+1];
    223   strcpy(buf, instore.str().data());
    224   return buf;
     219    ostringstream instore;
     220    instore << (*x);
     221    int n = strlen(instore.str().data());
     222    char* buf = new char[n+1];
     223    strcpy(buf, instore.str().data());
     224    return buf;
    225225}
    226226
    227227struct ZZX* ZZX_copy(struct ZZX* x) {
    228   return new ZZX(*x);
     228    return new ZZX(*x);
    229229}
    230230
    231231/* Sets ith coefficient of x to value.
    232232   AUTHOR: David Harvey (2006-06-08) */
    233233void ZZX_setitem_from_int(struct ZZX* x, long i, int value)
    234234{
    235   SetCoeff(*x, i, value);
     235    SetCoeff(*x, i, value);
    236236}
    237237
    238238/* Returns ith coefficient of x.
     
    240240   AUTHOR: David Harvey (2006-06-08) */
    241241int ZZX_getitem_as_int(struct ZZX* x, long i)
    242242{
    243   return ZZ_to_int(&coeff(*x, i));
     243    return ZZ_to_int(&coeff(*x, i));
    244244}
    245245
    246246/* Copies ith coefficient of x to output.
     
    254254
    255255struct ZZX* ZZX_div(struct ZZX* x, struct ZZX* y, int* divisible)
    256256{
    257   struct ZZX* z = new ZZX();
    258   *divisible = divide(*z, *x, *y);
    259   return z;
     257    struct ZZX* z = new ZZX();
     258    *divisible = divide(*z, *x, *y);
     259    return z;
    260260}
    261261
    262262
    263263
    264264void ZZX_quo_rem(struct ZZX* x, struct ZZX* other, struct ZZX** r, struct ZZX** q)
    265265{
    266   struct ZZX *qq = new ZZX(), *rr = new ZZX();
    267   DivRem(*qq, *rr, *x, *other);
    268   *r = rr; *q = qq;
     266    struct ZZX *qq = new ZZX(), *rr = new ZZX();
     267    DivRem(*qq, *rr, *x, *other);
     268    *r = rr; *q = qq;
    269269}
    270270
    271271
    272272struct ZZX* ZZX_square(struct ZZX* x)
    273273{
    274   struct ZZX* s = new ZZX();
    275   sqr(*s, *x);
    276   return s;
     274    struct ZZX* s = new ZZX();
     275    sqr(*s, *x);
     276    return s;
    277277}
    278278
    279279
    280280int ZZX_is_monic(struct ZZX* x)
    281281{
    282   return IsOne(LeadCoeff(*x));
     282    return IsOne(LeadCoeff(*x));
    283283}
    284284
    285285
    286286struct ZZX* ZZX_neg(struct ZZX* x)
    287287{
    288   struct ZZX* y = new ZZX();
    289   *y = -*x;
    290   return y;
     288    struct ZZX* y = new ZZX();
     289    *y = -*x;
     290    return y;
    291291}
    292292
    293293
    294294struct ZZX* ZZX_left_shift(struct ZZX* x, long n)
    295295{
    296   struct ZZX* y = new ZZX();
    297   LeftShift(*y, *x, n);
    298   return y;
     296    struct ZZX* y = new ZZX();
     297    LeftShift(*y, *x, n);
     298    return y;
    299299}
    300300
    301301
    302302struct ZZX* ZZX_right_shift(struct ZZX* x, long n)
    303303{
    304   struct ZZX* y = new ZZX();
    305   RightShift(*y, *x, n);
    306   return y;
     304    struct ZZX* y = new ZZX();
     305    RightShift(*y, *x, n);
     306    return y;
    307307}
    308308
    309309struct ZZX* ZZX_primitive_part(struct ZZX* x)
    310310{
    311   struct ZZX* p = new ZZX();
    312   PrimitivePart(*p, *x);
    313   return p;
     311    struct ZZX* p = new ZZX();
     312    PrimitivePart(*p, *x);
     313    return p;
    314314}
    315315
    316316
    317317void ZZX_pseudo_quo_rem(struct ZZX* x, struct ZZX* y, struct ZZX** r, struct ZZX** q)
    318318{
    319   *r = new ZZX();
    320   *q = new ZZX();
    321   PseudoDivRem(**q, **r, *x, *y);
     319    *r = new ZZX();
     320    *q = new ZZX();
     321    PseudoDivRem(**q, **r, *x, *y);
    322322}
    323323
    324324
    325325struct ZZX* ZZX_gcd(struct ZZX* x, struct ZZX* y)
    326326{
    327   struct ZZX* g = new ZZX();
    328   GCD(*g, *x, *y);
    329   return g;
     327    struct ZZX* g = new ZZX();
     328    GCD(*g, *x, *y);
     329    return g;
    330330}
    331331
    332332
    333 void ZZX_xgcd(struct ZZX* x, struct ZZX* y, struct ZZ** r, struct ZZX** s, \
    334               struct ZZX** t, int proof)
     333void ZZX_xgcd(struct ZZX* x, struct ZZX* y, struct ZZ** r, struct ZZX** s,
     334          struct ZZX** t, int proof)
    335335{
    336   *r = new ZZ();
    337   *s = new ZZX();
    338   *t = new ZZX();
    339   XGCD(**r, **s, **t, *x, *y, proof);
     336    *r = new ZZ();
     337    *s = new ZZX();
     338    *t = new ZZX();
     339    XGCD(**r, **s, **t, *x, *y, proof);
    340340}
    341341
    342342
    343343long ZZX_degree(struct ZZX* x)
    344344{
    345   return deg(*x);
     345    return deg(*x);
    346346}
    347347
    348348void ZZX_set_x(struct ZZX* x)
    349349{
    350   SetX(*x);
     350    SetX(*x);
    351351}
    352352
    353353
    354354int ZZX_is_x(struct ZZX* x)
    355355{
    356   return IsX(*x);
     356    return IsX(*x);
    357357}
    358358
    359359
    360360struct ZZX* ZZX_derivative(struct ZZX* x)
    361361{
    362   ZZX* d = new ZZX();
    363   diff(*d, *x);
    364   return d;
     362    ZZX* d = new ZZX();
     363    diff(*d, *x);
     364    return d;
    365365}
    366366
    367367
    368368struct ZZX* ZZX_reverse(struct ZZX* x)
    369369{
    370   ZZX* r = new ZZX();
    371   reverse(*r, *x);
    372   return r;
     370    ZZX* r = new ZZX();
     371    reverse(*r, *x);
     372    return r;
    373373}
    374374
    375375struct ZZX* ZZX_reverse_hi(struct ZZX* x, int hi)
    376376{
    377   ZZX* r = new ZZX();
    378   reverse(*r, *x, hi);
    379   return r;
     377    ZZX* r = new ZZX();
     378    reverse(*r, *x, hi);
     379    return r;
    380380}
    381381
    382382
    383383struct ZZX* ZZX_truncate(struct ZZX* x, long m)
    384384{
    385   ZZX* t = new ZZX();
    386   trunc(*t, *x, m);
    387   return t;
     385    ZZX* t = new ZZX();
     386    trunc(*t, *x, m);
     387    return t;
    388388}
    389389
    390390
    391391struct ZZX* ZZX_multiply_and_truncate(struct ZZX* x, struct ZZX* y, long m)
    392392{
    393   ZZX* t = new ZZX();
    394   MulTrunc(*t, *x, *y, m);
    395   return t;
     393    ZZX* t = new ZZX();
     394    MulTrunc(*t, *x, *y, m);
     395    return t;
    396396}
    397397
    398398
    399399struct ZZX* ZZX_square_and_truncate(struct ZZX* x, long m)
    400400{
    401   ZZX* t = new ZZX();
    402   SqrTrunc(*t, *x, m);
    403   return t;
     401    ZZX* t = new ZZX();
     402    SqrTrunc(*t, *x, m);
     403    return t;
    404404}
    405405
    406406
    407407struct ZZX* ZZX_invert_and_truncate(struct ZZX* x, long m)
    408408{
    409   ZZX* t = new ZZX();
    410   InvTrunc(*t, *x, m);
    411   return t;
     409    ZZX* t = new ZZX();
     410    InvTrunc(*t, *x, m);
     411    return t;
    412412}
    413413
    414414
    415415struct ZZX* ZZX_multiply_mod(struct ZZX* x, struct ZZX* y,  struct ZZX* modulus)
    416416{
    417   ZZX* p = new ZZX();
    418   MulMod(*p, *x, *y, *modulus);
    419   return p;
     417    ZZX* p = new ZZX();
     418    MulMod(*p, *x, *y, *modulus);
     419    return p;
    420420}
    421421
    422422
    423423struct ZZ* ZZX_trace_mod(struct ZZX* x, struct ZZX* y)
    424424{
    425   ZZ* p = new ZZ();
    426   TraceMod(*p, *x, *y);
    427   return p;
     425    ZZ* p = new ZZ();
     426    TraceMod(*p, *x, *y);
     427    return p;
    428428}
    429429
    430430
    431431char* ZZX_trace_list(struct ZZX* x)
    432432{
    433   vec_ZZ v;
    434   TraceVec(v, *x);
    435   ostringstream instore;
    436   instore << v;
    437   int n = strlen(instore.str().data());
    438   char* buf = new char[n+1];
    439   strcpy(buf, instore.str().data());
    440   return buf;
     433    vec_ZZ v;
     434    TraceVec(v, *x);
     435    ostringstream instore;
     436    instore << v;
     437    int n = strlen(instore.str().data());
     438    char* buf = new char[n+1];
     439    strcpy(buf, instore.str().data());
     440    return buf;
    441441}
    442442
    443443
    444444struct ZZ* ZZX_resultant(struct ZZX* x, struct ZZX* y, int proof)
    445445{
    446   ZZ* res = new ZZ();
    447   resultant(*res, *x, *y, proof);
    448   return res;
     446    ZZ* res = new ZZ();
     447    resultant(*res, *x, *y, proof);
     448    return res;
    449449}
    450450
    451451
    452452struct ZZ* ZZX_norm_mod(struct ZZX* x, struct ZZX* y, int proof)
    453453{
    454   ZZ* res = new ZZ();
    455   NormMod(*res, *x, *y, proof);
    456   return res;
     454    ZZ* res = new ZZ();
     455    NormMod(*res, *x, *y, proof);
     456    return res;
    457457}
    458458
    459459
    460460struct ZZ* ZZX_discriminant(struct ZZX* x, int proof)
    461461{
    462   ZZ* d = new ZZ();
    463   discriminant(*d, *x, proof);
    464   return d;
     462    ZZ* d = new ZZ();
     463    discriminant(*d, *x, proof);
     464    return d;
    465465}
    466466
    467467
    468468struct ZZX* ZZX_charpoly_mod(struct ZZX* x, struct ZZX* y, int proof)
    469469{
    470   ZZX* f = new ZZX();
    471   CharPolyMod(*f, *x, *y, proof);
    472   return f;
     470    ZZX* f = new ZZX();
     471    CharPolyMod(*f, *x, *y, proof);
     472    return f;
    473473}
    474474
    475475
    476476struct ZZX* ZZX_minpoly_mod(struct ZZX* x, struct ZZX* y)
    477477{
    478   ZZX* f = new ZZX();
    479   MinPolyMod(*f, *x, *y);
    480   return f;
     478    ZZX* f = new ZZX();
     479    MinPolyMod(*f, *x, *y);
     480    return f;
    481481}
    482482
    483483
    484484void ZZX_clear(struct ZZX* x)
    485485{
    486   clear(*x);
     486    clear(*x);
    487487}
    488488
    489489
    490490void ZZX_preallocate_space(struct ZZX* x, long n)
    491491{
    492   x->SetMaxLength(n);
     492    x->SetMaxLength(n);
    493493}
    494494
    495495/*
    496496EXTERN struct ZZ* ZZX_polyeval(struct ZZX* f, struct ZZ* a)
    497497{
    498   ZZ* b = new ZZ();
    499   *b = PolyEval(*f, *a);
    500   return b;
     498    ZZ* b = new ZZ();
     499    *b = PolyEval(*f, *a);
     500    return b;
    501501}
    502502*/
    503503
    504504void ZZX_squarefree_decomposition(struct ZZX*** v, long** e, long* n, struct ZZX* x)
    505505{
    506   vec_pair_ZZX_long factors;
    507   SquareFreeDecomp(factors, *x);
    508   *n = factors.length();
    509   *v = (ZZX**) malloc(sizeof(ZZX*) * (*n));
    510   *e = (long*) malloc(sizeof(long) * (*n));
    511   for (long i = 0; i < (*n); i++) {
    512     (*v)[i] = new ZZX(factors[i].a);
    513     (*e)[i] = factors[i].b;
    514   }
     506    vec_pair_ZZX_long factors;
     507    SquareFreeDecomp(factors, *x);
     508    *n = factors.length();
     509    *v = (ZZX**) malloc(sizeof(ZZX*) * (*n));
     510    *e = (long*) malloc(sizeof(long) * (*n));
     511    for (long i = 0; i < (*n); i++) {
     512        (*v)[i] = new ZZX(factors[i].a);
     513        (*e)[i] = factors[i].b;
     514    }
    515515}
    516516
    517517///////////////////////////////////////////////
     
    753753
    754754char* ZZ_pX_trace_list(struct ZZ_pX* x)
    755755{
    756   vec_ZZ_p v;
    757   TraceVec(v, *x);
    758   ostringstream instore;
    759   instore << v;
    760   int n = strlen(instore.str().data());
    761   char* buf = new char[n+1];
    762   strcpy(buf, instore.str().data());
    763   return buf;
     756    vec_ZZ_p v;
     757    TraceVec(v, *x);
     758    ostringstream instore;
     759    instore << v;
     760    int n = strlen(instore.str().data());
     761    char* buf = new char[n+1];
     762    strcpy(buf, instore.str().data());
     763    return buf;
    764764}
    765765
    766766
     
    810810
    811811void ZZ_pX_factor(struct ZZ_pX*** v, long** e, long* n, struct ZZ_pX* x, long verbose)
    812812{
    813   long i;
    814   vec_pair_ZZ_pX_long factors;
    815   berlekamp(factors, *x, verbose);
    816   *n = factors.length();
    817   *v = (ZZ_pX**) malloc(sizeof(ZZ_pX*) * (*n));
    818   *e = (long*) malloc(sizeof(long)*(*n));
    819   for (i=0; i<(*n); i++) {
    820     (*v)[i] = new ZZ_pX(factors[i].a);
    821     (*e)[i] = factors[i].b;
    822   }
     813    long i;
     814    vec_pair_ZZ_pX_long factors;
     815    berlekamp(factors, *x, verbose);
     816    *n = factors.length();
     817    *v = (ZZ_pX**) malloc(sizeof(ZZ_pX*) * (*n));
     818    *e = (long*) malloc(sizeof(long)*(*n));
     819    for (i=0; i<(*n); i++) {
     820        (*v)[i] = new ZZ_pX(factors[i].a);
     821        (*e)[i] = factors[i].b;
     822    }
    823823}
    824824
    825825void ZZ_pX_linear_roots(struct ZZ_p*** v, long* n, struct ZZ_pX* f)
    826826{
    827   long i;
    828   // printf("1\n");
    829   vec_ZZ_p w;
    830   FindRoots(w, *f);
    831   // printf("2\n");
    832   *n = w.length();
    833   //   printf("3 %d\n",*n);
    834   (*v) = (ZZ_p**) malloc(sizeof(ZZ_p*)*(*n));
    835   for (i=0; i<(*n); i++) {
    836     (*v)[i] = new ZZ_p(w[i]);
    837   }
     827    long i;
     828    vec_ZZ_p w;
     829    FindRoots(w, *f);
     830    *n = w.length();
     831    (*v) = (ZZ_p**) malloc(sizeof(ZZ_p*)*(*n));
     832    for (i=0; i<(*n); i++) {
     833        (*v)[i] = new ZZ_p(w[i]);
     834    }
    838835}
    839836
    840837/////////// ZZ_pE //////////////
    841838
    842839struct ZZ_pX ZZ_pE_to_ZZ_pX(struct ZZ_pE x)
    843840{
    844   return ZZ_pX(rep(x));
     841    return ZZ_pX(rep(x));
    845842}
    846843
    847844
     
    849846//////// mat_ZZ //////////
    850847
    851848void mat_ZZ_SetDims(struct mat_ZZ* mZZ, long nrows, long ncols){
    852   mZZ->SetDims(nrows, ncols);
     849    mZZ->SetDims(nrows, ncols);
    853850}
    854851
    855852struct mat_ZZ* mat_ZZ_pow(const struct mat_ZZ* x, long e)
    856853{
    857   mat_ZZ *z = new mat_ZZ();
    858   power(*z, *x, e);
    859   return z;
     854    mat_ZZ *z = new mat_ZZ();
     855    power(*z, *x, e);
     856    return z;
    860857}
    861858
    862859long mat_ZZ_nrows(const struct mat_ZZ* x)
    863860{
    864   return x->NumRows();
     861    return x->NumRows();
    865862}
    866863
    867864
    868865long mat_ZZ_ncols(const struct mat_ZZ* x)
    869866{
    870   return x->NumCols();
     867    return x->NumCols();
    871868}
    872869
    873870void mat_ZZ_setitem(struct mat_ZZ* x, int i, int j, const struct ZZ* z)
    874871{
    875   (*x)[i][j] = *z;
     872    (*x)[i][j] = *z;
    876873 
    877874}
    878875
    879876struct ZZ* mat_ZZ_getitem(const struct mat_ZZ* x, int i, int j)
    880877{
    881   return new ZZ((*x)(i,j));
     878    return new ZZ((*x)(i,j));
    882879}
    883880
    884881struct ZZ* mat_ZZ_determinant(const struct mat_ZZ* x, long deterministic)
    885882{
    886   ZZ* d = new ZZ();
    887   determinant(*d, *x, deterministic);
    888   return d;
     883    ZZ* d = new ZZ();
     884    determinant(*d, *x, deterministic);
     885    return d;
    889886}
    890887
    891888struct mat_ZZ* mat_ZZ_HNF(const struct mat_ZZ* A, const struct ZZ* D)
    892889{
    893   struct mat_ZZ* W = new mat_ZZ();
    894   HNF(*W, *A, *D);
    895   return W;
     890    struct mat_ZZ* W = new mat_ZZ();
     891    HNF(*W, *A, *D);
     892    return W;
    896893}
    897894
    898895long mat_ZZ_LLL(struct ZZ **det, struct mat_ZZ *x, long a, long b, long verbose)
    899896{
    900   *det = new ZZ();
    901   return LLL(**det,*x,a,b,verbose);
     897    *det = new ZZ();
     898    return LLL(**det,*x,a,b,verbose);
    902899}
    903900
    904901long mat_ZZ_LLL_U(struct ZZ **det, struct mat_ZZ *x, struct mat_ZZ *U, long a, long b, long verbose)
    905902{
    906   *det = new ZZ();
    907   return LLL(**det,*x,*U,a,b,verbose);
     903    *det = new ZZ();
     904    return LLL(**det,*x,*U,a,b,verbose);
    908905}
    909906
    910907
    911908struct ZZX* mat_ZZ_charpoly(const struct mat_ZZ* A)
    912909{
    913   ZZX* f = new ZZX();
    914   CharPoly(*f, *A);
    915   return f;
     910    ZZX* f = new ZZX();
     911    CharPoly(*f, *A);
     912    return f;
    916913}
    917914
    918915/**
     
    921918
    922919GF2EContext* GF2EContext_construct(void *mem, const GF2X *p)
    923920{
    924   return new(mem) GF2EContext(*p);
     921    return new(mem) GF2EContext(*p);
    925922}
    926923
    927924
    928925GF2EContext* GF2EContext_new(const GF2X *p)
    929926{
    930   return new GF2EContext(*p);
     927    return new GF2EContext(*p);
    931928}
    932929
    933930
    934931void mat_GF2E_setitem(struct mat_GF2E* x, int i, int j, const struct GF2E* z)
    935932{
    936   (*x)[i][j] = *z;
     933    (*x)[i][j] = *z;
    937934}
    938935
    939936void mat_GF2_setitem(struct mat_GF2* x, int i, int j, const struct GF2* z)
    940937{
    941   (*x)[i][j] = *z;
     938    (*x)[i][j] = *z;
    942939}
    943940
    944941/**
     
    947944
    948945ZZ_pContext* ZZ_pContext_new(ZZ *p)
    949946{
    950         return new ZZ_pContext(*p);
     947    return new ZZ_pContext(*p);
    951948}
    952949
    953950ZZ_pContext* ZZ_pContext_construct(void *mem, ZZ *p)
    954951{
    955         return new(mem) ZZ_pContext(*p);
     952    return new(mem) ZZ_pContext(*p);
    956953}
    957954
    958955void ZZ_pContext_restore(ZZ_pContext *ctx)
    959956{
    960         ctx->restore();
     957    ctx->restore();
    961958}
    962959
    963960// Functions for using ZZ_pX's for p-adic extensions
     
    11961193
    11971194ZZ_pEContext* ZZ_pEContext_new(ZZ_pX *f)
    11981195{
    1199         return new ZZ_pEContext(*f);
     1196    return new ZZ_pEContext(*f);
    12001197}
    12011198
    12021199ZZ_pEContext* ZZ_pEContext_construct(void *mem, ZZ_pX *f)
    12031200{
    1204         return new(mem) ZZ_pEContext(*f);
     1201    return new(mem) ZZ_pEContext(*f);
    12051202}
    12061203
    12071204void ZZ_pEContext_restore(ZZ_pEContext *ctx)
    12081205{
    1209         ctx->restore();
     1206    ctx->restore();
    12101207}
    12111208
    12121209/**
     
    12151212
    12161213zz_pContext* zz_pContext_new(long p)
    12171214{
    1218         return new zz_pContext(p);
     1215    return new zz_pContext(p);
    12191216}
    12201217
    12211218zz_pContext* zz_pContext_construct(void *mem, long p)
    12221219{
    1223         return new(mem) zz_pContext(p);
     1220    return new(mem) zz_pContext(p);
    12241221}
    12251222
    12261223void zz_pContext_restore(zz_pContext *ctx)
    12271224{
    1228         ctx->restore();
     1225    ctx->restore();
    12291226}
  • doc/common/themes/sage/layout.html

    diff --git a/doc/common/themes/sage/layout.html b/doc/common/themes/sage/layout.html
    a b  
    160160            <h3>{{ _('Quick search') }}</h3>
    161161              <form class="search" action="{{ pathto('search') }}" method="get">
    162162                <input type="text" name="q" size="18" />
    163                 <!-- The shading of the "Go" button should be consistent -->
    164                 <!-- with the colour of the header and footer. See the file -->
    165                 <!-- doc/common/themes/sage/theme.conf for colours used by -->
    166                 <!-- the Sage theme. -->
     163                <!-- The shading of the "Go" button should be consistent -->
     164                <!-- with the colour of the header and footer. See the file -->
     165                <!-- doc/common/themes/sage/theme.conf for colours used by -->
     166                <!-- the Sage theme. -->
    167167                <input type="submit" style="background-color: #B8B9F6" value="{{ _('Go') }}" />
    168168                <input type="hidden" name="check_keywords" value="yes" />
    169169                <input type="hidden" name="area" value="default" />
  • doc/en/developer/coding_in_python.rst

    diff --git a/doc/en/developer/coding_in_python.rst b/doc/en/developer/coding_in_python.rst
    a b  
    139139            sage: float(pi)
    140140            3.1415926535897931
    141141        """
    142         ...
     142        ...
    143143        def _repr_(self):
    144144            return "pi"
    145145
     
    496496
    497497    class Foo(SageObject):
    498498        def terrible_idea(self):
    499             return 1
     499            return 1
    500500        def bad_name(self):
    501             return 1
     501            return 1
    502502        def f(self, weird_keyword=True):
    503             return 1
     503            return 1
    504504
    505505You note that the ``terrible_idea()`` method does not make any sense,
    506506and should be removed altogether. You open the trac ticket number 3333
    507507(say), and replace the code with::
    508508
    509509        def terrible_idea(self):
    510             from sage.misc.superseded import deprecation
    511             deprecation(3333, 'You can just call f() instead')
    512             return 1
     510            from sage.misc.superseded import deprecation
     511            deprecation(3333, 'You can just call f() instead')
     512            return 1
    513513
    514514Later, you come up with a much better name for the second method. You
    515515open the trac ticket number 4444, and replace it with::
    516516
    517517        def much_better_name(self):
    518             return 1
     518            return 1
    519519
    520         bad_name = deprecated_function_alias(4444, much_better_name)
     520        bad_name = deprecated_function_alias(4444, much_better_name)
    521521
    522522Finally, you like the ``f()`` method name but you don't like the
    523523``weird_keyword`` name. You fix this by opening the trac ticket 5555,
     
    525525
    526526        @rename_keyword(deprecation=5555, weird_keyword='nice_keyword')
    527527        def f(self, nice_keyword=True):
    528             return 1
     528            return 1
    529529
    530530Now, any user that still relies on the deprecated functionality will
    531531be informed that this is about to change, yet the deprecated commands
     
    538538    sage: class Foo(SageObject):
    539539    ...
    540540    ...     def terrible_idea(self):
    541     ...             deprecation(3333, 'You can just call f() instead')
    542     ...             return 1
     541    ...         deprecation(3333, 'You can just call f() instead')
     542    ...         return 1
    543543    ...
    544544    ...     def much_better_name(self):
    545     ...             return 1
     545    ...         return 1
    546546    ...
    547547    ...     bad_name = deprecated_function_alias(4444, much_better_name)
    548548    ...
    549549    ...     @rename_keyword(deprecation=5555, weird_keyword='nice_keyword')
    550550    ...     def f(self, nice_keyword=True):
    551     ...             return 1
     551    ...         return 1
    552552
    553553    sage: foo = Foo()
    554554    sage: foo.terrible_idea()
     
    565565    doctest:...: DeprecationWarning: use the option 'nice_keyword' instead of 'weird_keyword'
    566566    See http://trac.sagemath.org/5555 for details.
    567567    1
    568  
    569568
    570569Editing existing files
    571570======================
  • doc/en/developer/conventions.rst

    diff --git a/doc/en/developer/conventions.rst b/doc/en/developer/conventions.rst
    a b  
    418418      .. TODO::
    419419
    420420          Improve further function ``have_fresh_beers`` using algorithm
    421           ``buy_a_better_fridge``::
     421          ``buy_a_better_fridge``::
    422422
    423               sage: have_fresh_beers('Bière de l'Yvette') # todo: not implemented
    424               Enjoy !
     423              sage: have_fresh_beers('Bière de l'Yvette') # todo: not implemented
     424              Enjoy !
    425425
    426426- A REFERENCES block to list books or papers (optional). This block serves
    427427  a similar purpose to a list of references in a research paper, or a
     
    468468
    469469        The point as a tuple.
    470470
    471             .. SEEALSO::
     471        .. SEEALSO::
    472472
    473473            :func:`line`
    474474
     
    527527        If it contains multiple lines, all the lines after the
    528528        first need to be indented.
    529529
    530         :type x: integer; default 1
     530        :type x: integer; default 1
    531531
    532         :param y: the ...
     532        :param y: the ...
    533533
    534         :type y: integer; default 2
     534        :type y: integer; default 2
    535535
    536         :returns: the ...
     536        :returns: the ...
    537537
    538         :rtype: integer, the return type
     538        :rtype: integer, the return type
    539539
    540540You are strongly encouraged to:
    541541
     
    624624
    625625.. WARNING::
    626626
    627    Functions whose names start with an underscore do not currently
    628    appear in the reference manual, so avoid putting crucial
    629    documentation in their docstrings. In particular, if you are
    630    defining a class, you might put a long informative docstring after
    631    the class definition, not for the ``__init__`` method. For example,
    632    from the file ``SAGE_ROOT/devel/sage/sage/crypto/classical.py``:
     627    Functions whose names start with an underscore do not currently
     628    appear in the reference manual, so avoid putting crucial
     629    documentation in their docstrings. In particular, if you are
     630    defining a class, you might put a long informative docstring after
     631    the class definition, not for the ``__init__`` method. For example,
     632    from the file ``SAGE_ROOT/devel/sage/sage/crypto/classical.py``::
    633633
    634    ::
     634        class HillCryptosystem(SymmetricKeyCryptosystem):
     635            """
     636            Create a Hill cryptosystem defined by the `m` x `m` matrix space
     637            over `\mathbf{Z} / N \mathbf{Z}`, where `N` is the alphabet size of
     638            the string monoid ``S``.
    635639
    636     class HillCryptosystem(SymmetricKeyCryptosystem):
     640            INPUT:
     641
     642            - ``S`` -- a string monoid over some alphabet
     643
     644            - ``m`` -- a positive integer; the block length of matrices that
     645              specify block permutations
     646
     647            OUTPUT:
     648
     649            - A Hill cryptosystem of block length ``m`` over the alphabet ``S``.
     650
     651            EXAMPLES::
     652
     653                sage: S = AlphabeticStrings()
     654                sage: E = HillCryptosystem(S,3)
     655                sage: E
     656                Hill cryptosystem on Free alphabetic string monoid on A-Z of block length 3
    637657        """
    638         Create a Hill cryptosystem defined by the `m` x `m` matrix space
    639         over `\mathbf{Z} / N \mathbf{Z}`, where `N` is the alphabet size of
    640         the string monoid ``S``.
    641658
    642         INPUT:
    643 
    644         - ``S`` -- a string monoid over some alphabet
    645 
    646         - ``m`` -- a positive integer; the block length of matrices that
    647           specify block permutations
    648 
    649         OUTPUT:
    650 
    651         - A Hill cryptosystem of block length ``m`` over the alphabet ``S``.
    652 
    653         EXAMPLES::
    654 
    655             sage: S = AlphabeticStrings()
    656             sage: E = HillCryptosystem(S,3)
    657             sage: E
    658             Hill cryptosystem on Free alphabetic string monoid on A-Z of block length 3
    659         """
    660 
    661    and so on, while the ``__init__`` method starts like this::
     659    and so on, while the ``__init__`` method starts like this::
    662660
    663661        def __init__(self, S, m):
    664662            """
    665663            See ``HillCryptosystem`` for full documentation.
    666664
    667             EXAMPLES::
    668             ...
    669             """
     665        EXAMPLES::
     666        ...
     667        """
    670668
    671    Note also that the first docstring is printed if users type
    672    "HillCryptosystem?" at the "sage:" prompt.
     669    Note also that the first docstring is printed if users type
     670    "HillCryptosystem?" at the "sage:" prompt.
    673671
    674    (Before Sage 3.4, the reference manual used to include methods
    675    starting with underscores, so you will probably find many examples
    676    in the code which don't follow this advice...)
     672    (Before Sage 3.4, the reference manual used to include methods
     673    starting with underscores, so you will probably find many examples
     674    in the code which don't follow this advice...)
    677675
    678676
    679677Automatic testing
     
    892890     The following should yield 4.  See :trac:`2`. ::
    893891
    894892        sage: 2+2 # optional: bug
    895         5
     893        5
    896894
    897895  Then the doctest will be skipped by default, but could be revealed
    898896  by running ``sage -t --only-optional=bug ...``.  (A doctest marked
  • doc/en/thematic_tutorials/lie/weyl_character_ring.rst

    diff --git a/doc/en/thematic_tutorials/lie/weyl_character_ring.rst b/doc/en/thematic_tutorials/lie/weyl_character_ring.rst
    a b  
    552552
    553553::
    554554
    555         sage: A2=WeylCharacterRing("A2",style="coroots")
    556         sage: ad = A2(1,1)
    557         sage: [ad.symmetric_power(k).invariant_degree() for k in [0..6]]
    558         [1, 0, 1, 1, 1, 1, 2]
    559         sage: [ad.exterior_power(k).invariant_degree() for k in [0..6]]
    560         [1, 0, 0, 1, 0, 1, 0]
     555    sage: A2=WeylCharacterRing("A2",style="coroots")
     556    sage: ad = A2(1,1)
     557    sage: [ad.symmetric_power(k).invariant_degree() for k in [0..6]]
     558    [1, 0, 1, 1, 1, 1, 2]
     559    sage: [ad.exterior_power(k).invariant_degree() for k in [0..6]]
     560    [1, 0, 0, 1, 0, 1, 0]
    561561
    562562If we want the multiplicity of some other representation, we may
    563563obtain that using the method ``multiplicity``::
     
    566566
    567567::
    568568
    569         sage: (ad^3).multiplicity(ad)
    570         8
     569    sage: (ad^3).multiplicity(ad)
     570    8
    571571
  • doc/en/website/templates/index.html

    diff --git a/doc/en/website/templates/index.html b/doc/en/website/templates/index.html
    a b  
    8888
    8989    <td width="50%">
    9090      <p class="biglink">
    91         <a class="biglink" href="thematic_tutorials/index.html">
    92           Thematic Tutorials
    93         </a>
    94         <a title="Download PDF" class="pdf"
    95            href="../../pdf/en/thematic_tutorials/thematic_tutorials.pdf">
    96           <img class="icon" src="_static/pdf.png"></img>
    97         </a>
    98         <br>
    99         <span class="linkdescr">
    100           A collection of in-depth tutorials on specific topics. These
    101           thematic tutorials are designed to help you get started on
    102           Sage functionalities relating to topics such as coding
    103           theory, combinatorics, cryptography, functional programming,
    104           group theory, linear programming, etc. If you feel
    105           uncomfortable consulting the reference manual, you are
    106           encouraged to browse through these thematic tutorials.
    107         </span>
     91        <a class="biglink" href="thematic_tutorials/index.html">
     92          Thematic Tutorials
     93        </a>
     94        <a title="Download PDF" class="pdf"
     95           href="../../pdf/en/thematic_tutorials/thematic_tutorials.pdf">
     96          <img class="icon" src="_static/pdf.png"></img>
     97        </a>
     98        <br>
     99        <span class="linkdescr">
     100          A collection of in-depth tutorials on specific topics. These
     101          thematic tutorials are designed to help you get started on
     102          Sage functionalities relating to topics such as coding
     103          theory, combinatorics, cryptography, functional programming,
     104          group theory, linear programming, etc. If you feel
     105          uncomfortable consulting the reference manual, you are
     106          encouraged to browse through these thematic tutorials.
     107        </span>
    108108      </p>
    109109    </td>
    110110
  • sage/combinat/matrices/dancing_links_c.h

    diff --git a/sage/combinat/matrices/dancing_links_c.h b/sage/combinat/matrices/dancing_links_c.h
    a b  
    3535template <class T>
    3636string to_string(T x)
    3737{
    38         std::ostringstream stream_out;
    39         stream_out << x;
    40         return stream_out.str();
     38    std::ostringstream stream_out;
     39    stream_out << x;
     40    return stream_out.str();
    4141}
    4242
    4343typedef struct node_struct {
  • sage/combinat/partitions_c.cc

    diff --git a/sage/combinat/partitions_c.cc b/sage/combinat/partitions_c.cc
    a b  
    2424 *      a(n,k) = sum_{h=1, (h,k) = 1}^k exp(\pi i s(h,k) - 2 \pi i h  n / k)
    2525 *
    2626 *      and
    27  *     
     27 *
    2828 *      f_n(k) = \pi sqrt{2} cosh(A_n/(sqrt{3}*k))/(B_n*k) - sinh(C_n/k)/D_n;
    2929 *
    3030 *      where
     
    4646 *          -- Search source code for other TODO comments.
    4747 *
    4848 *      OTHER CREDITS:
    49  *     
     49 *
    5050 *      I looked source code written by Ralf Stephan, currently available at
    5151 *
    5252 *              http://www.ark.in-berlin.de/part.c
     
    8080 *      GNU General Public License for more details.
    8181 *
    8282 *      You should have received a copy of the GNU General Public License
    83  *      along with this program.  If not, see <http://www.gnu.org/licenses/>. 
     83 *      along with this program.  If not, see <http://www.gnu.org/licenses/>.
    8484 */
    8585
    8686
     
    126126
    127127const bool debug_precision = false;                             // If true, output information that might be useful for
    128128                                                                // debugging the precision setting code.
    129                                        
     129
    130130const bool debugs = false;                                      // If true, output informaiton that might be useful for
    131131                                                                // debugging the various s() functions.
    132132const bool debugf = false;                                      // Same for the f() functions.
     
    155155                                                                            // The assumed precision of a long double.
    156156                                                                            // Note: On many systems double_precision = long_double_precision. This is OK, as
    157157                                                                            // the long double stage of the computation will just be skipped.
    158                                                                             // 
     158                                                                            //
    159159                                                                            // NOTE: If long_double_precision > dd_precision, then, again, long doubles
    160160                                                                            //      will not ever be used. It would be nice if this were fixed.
    161161
     
    195195                                                                                        // because we don't know how much precision we will need
    196196                                                                                        // until we know for what n we are computing p(n).
    197197
    198 mpfr_t mp_A, mp_B, mp_C, mp_D;                                                          // These "constants" all depend on n, and 
     198mpfr_t mp_A, mp_B, mp_C, mp_D;                                                          // These "constants" all depend on n, and
    199199double d_A, d_B, d_C, d_D;                                                              //
    200200long double ld_A, ld_B, ld_C, ld_D;                                                     //
    201201
     
    237237 *
    238238 ****************************************************************************/
    239239
    240 // 
     240//
    241241// A listing of the main functions, in "conceptual" order.
    242242//
    243243
     
    331331                                                                                // have been verified by other programs, and
    332332                                                                                // also on known congruences for p(n)
    333333
    334 static void cospi (mpfr_t res, mpfr_t x);                                       // Puts cos(pi x) in result. This is not currently   
     334static void cospi (mpfr_t res, mpfr_t x);                                       // Puts cos(pi x) in result. This is not currently
    335335                                                                                // used, but should be faster than using mpfr_cos,
    336336                                                                                // according to comments in the file part.c
    337                                                                                 // mentioned in the introduction.   
     337                                                                                // mentioned in the introduction.
    338338                                                                                // test
    339339
    340340static int grab_last_digits(char * output, int n, mpfr_t x);                    // Might be useful for debugging, but
     
    355355 *
    356356 **********************************************************************/
    357357
    358 // The following function can be useful for debugging in come circumstances, but should not be used for anything else 
     358// The following function can be useful for debugging in come circumstances, but should not be used for anything else
    359359// unless it is rewritten.
    360360int grab_last_digits(char * output, int n, mpfr_t x) {
    361361    // fill output with the n digits of x that occur
     
    373373    int retval;
    374374
    375375    if(e > 0) {
    376         strncpy(output, temp + e - n, n);   
     376        strncpy(output, temp + e - n, n);
    377377        retval =  strlen(temp + e);
    378378    }
    379379    else {
     
    393393unsigned int compute_initial_precision(unsigned int n) {
    394394    // We just want to know how many bits we will need to
    395395    // compute to get an accurate answer.
    396    
     396
    397397    // We know that
    398398    //
    399399    //          p(n) ~ exp(pi * sqrt(2n/3))/(4n sqrt(3)),
     
    405405    // that the TOTAL ERROR in all computations is < (something small).
    406406    // This needs to be worked out carefully. EXTRA = log(n)/log(2) + 3
    407407    // is probably good enough, and is convenient...
    408     // 
     408    //
    409409    // but we really need:
    410410    //
    411     //                  p(n) < something 
     411    //                  p(n) < something
    412412    //
    413413    // to be sure that we compute the correct answer
    414414
    415415    unsigned int result = (unsigned int)(ceil(3.1415926535897931 * sqrt(2.0 * double(n)/ 3.0) / log(2))) + 3;
    416416    if(debug) cout << "Using initial precision of " << result << " bits." << endl;
    417    
     417
    418418    if(result > min_precision) {
    419419        return result;
    420420    }
    421    
     421
    422422    else return min_precision;
    423423
    424424}
    425425
    426426unsigned int compute_current_precision(unsigned int n, unsigned int N, unsigned int extra = 0) {
    427427    // Roughly, we compute
    428     // 
     428    //
    429429    //      log(A/sqrt(N) + B*sqrt(N/(n-1))*sinh(C * sqrt(n) / N) / log(2)
    430430    //
    431431    // where A, B, and C are the constants listed below. These error bounds
    432432    // are given in the paper by Rademacher listed at the top of this file.
    433433    // We then return this + extra, if extra != 0. If extra == 0, return with
    434434    // what is probably way more extra precision than is needed.
    435     // 
     435    //
    436436    // extra should probably have been set by a call to compute_extra_precision()
    437437    // before this function was called.
    438438
     
    448448    mpfr_init2(A, 32);
    449449    mpfr_init2(B, 32);
    450450    mpfr_init2(C, 32);
    451    
     451
    452452    mpfr_set_d(A,1.11431833485164,round_mode);
    453453    mpfr_set_d(B,0.059238439175445,round_mode);
    454454    mpfr_set_d(C,2.5650996603238,round_mode);
     
    472472    mpfr_set_ui(t2, N, round_mode);                       // t2 = N
    473473    mpfr_div_ui(t2, t2, n-1, round_mode);                 // t2 = N/(n-1)
    474474    mpfr_sqrt(t2, t2, round_mode);                        // t2 = sqrt( ditto )
    475    
     475
    476476    mpfr_mul(t1, t1, t2, round_mode);                     // t1 = B * sqrt(N/(n-1)) * sinh(C * sqrt(n)/N)
    477477
    478478    mpfr_add(error, error, t1, round_mode);               // error = (ERROR ESTIMATE)
    479                                                          
     479
    480480    unsigned int p = mpfr_get_exp(error);                 // I am not 100% certain that this always does the right thing.
    481481                                                          // (It should be the case that p is now the number of bits
    482482                                                          // required to hold the integer part of the error.)
    483    
     483
    484484
    485485    if(extra == 0) {
    486486        p = p + (unsigned int)ceil(log(n)/log(2));        // This is a stupid case to fall back on,
     
    529529    // in absolute value, and then return the extra precision
    530530    // that will guarantee that the accumulated error after computing
    531531    // that number of steps will be less than .5 - error.
    532     // 
     532    //
    533533    // How this works:
    534534    // We first need to figure out how many terms of the series we are going to
    535535    // need to compute. That is, we need to know how large k needs to be
     
    538538    // compute_current_precision() until we know that the error will be
    539539    // small enough to call compute_remainder(). Then we just call compute_remainder()
    540540    // until the error is small enough.
    541     // 
     541    //
    542542    // Now that we know how many terms we will need to compute, k, we compute
    543543    // the number of bits required to accurately store (.5 - error)/k. This ensures
    544544    // that the total error introduced when we add up all k terms of the sum
     
    559559                                                          // we end up doing a bunch of arithmetic operations, and if
    560560                                                          // we want the result of those operations to be accurate
    561561                                                          // within (.5 - error)/k, then we need that function to use
    562                                                           // a slightly higher working precision, which should be 
     562                                                          // a slightly higher working precision, which should be
    563563                                                          // independent of n.
    564564                                                          // TODO:
    565565                                                          // Extensive trial and error has found 3 to be the smallest value
     
    567567                                                          // be safe, we use 5 extra bits.
    568568                                                          // (Extensive trial and error means compiling this file to get
    569569                                                          // a.out and then running './a.out testforever' for a few hours.)
    570                                                          
    571                                                          
    572                                                        
     570
     571
     572
    573573
    574574    return bits;
    575575}
     
    650650    // NOTE: Calls to this function must be paired with calls to clear_constants()
    651651    static bool init = false;
    652652    mp_prec_t p = prec;
    653    
     653
    654654    mpfr_init2(mp_one_over_12,p); mpfr_init2(mp_one_over_24,p); mpfr_init2(mp_sqrt2,p); mpfr_init2(mp_sqrt3,p); mpfr_init2(mp_pi,p);
    655655    mpfr_init2(mp_A,p); mpfr_init2(mp_B,p); mpfr_init2(mp_C,p); mpfr_init2(mp_D,p); mpfr_init2(fourth, p); mpfr_init2(half, p);
    656    
     656
    657657    init = true;
    658    
     658
    659659    mpfr_set_ui(mp_one_over_12, 1, round_mode);                             // mp_one_over_12 = 1/12
    660660    mpfr_div_ui(mp_one_over_12, mp_one_over_12, 12, round_mode);            //
    661661
     
    701701    mpfr_div(mp_C, mp_C, mp_sqrt3, round_mode);                             // mp_C = sqrt(2) * pi * sqrt(n - 1/24) / sqrt3
    702702    //------------------------------------------------------------------------
    703703
    704    
     704
    705705    //mp_D = 2.0 * (n - 1.0/24.0) * sqrt(n - 1.0/24.0);-----------------------
    706706    mpfr_set_ui(mp_D, 2, round_mode);                                       // mp_D = 2
    707707    mpfr_mul(mp_D, mp_D, n_minus, round_mode);                              // mp_D = 2 * (n - 1/24)
     
    732732}
    733733
    734734void initialize_mpz_and_mpq_variables() {
    735     /* 
     735    /*
    736736     * We use a few mpz_t and mpq_t variables which need to be initialized
    737737     * before they can be used. Initialization and clearing take some
    738738     * time, so we initialize just once in this function, and clear in another.
     
    758758    // notice that this doesn't use n - the "constants"
    759759    // A, B, C, and D depend on n, but they are precomputed
    760760    // once before this function is called.
    761    
     761
    762762    //result =  pi * sqrt(2) * cosh(A/(sqrt(3)*k))/(B*k) - sinh(C/k)/D;
    763763
    764764                                                                    //
    765765    mpfr_set(result, mp_pi, round_mode);                            // result = pi
    766    
     766
    767767    mpfr_mul(result, result, mp_sqrt2, round_mode);                 // result = sqrt(2) * pi
    768768                                                                    //
    769769    mpfr_div(tempf1, mp_A, mp_sqrt3, round_mode);                   // temp1 = mp_A/sqrt(3)
     
    793793//       part of s(h,k)/2. It may be possible to make use of this somehow.
    794794//
    795795void q_s(mpq_t result, unsigned int h, unsigned int k) {
    796    
     796
    797797    if(k < 3) {
    798798        mpq_set_ui(result, 0, 1);
    799799        return;
     
    829829    //    return;
    830830    //}
    831831
    832 /*   
     832/*
    833833
    834834    // if k % h == 1, then
    835835    //
     
    858858    //if(k % h == 2) {
    859859    //}
    860860*/
    861    
    862861
    863862
    864    
     863
     864
    865865    mpq_set_ui(result, 0, 1);                             // result = 0
    866866
    867867    int r1 = k;
     
    878878            mpq_set_ui(qtemps, r1 * r1 + r2 * r2 + 1, r1 * r2);
    879879        }
    880880        //mpq_canonicalize(qtemps);
    881        
     881
    882882        if(n % 2 == 0){                                             //
    883883            mpq_add(result, result, qtemps);                        // result += temp1;
    884884        }                                                           //
     
    893893
    894894    mpq_set_ui(qtemps, 1, 12);
    895895    mpq_mul(result, result, qtemps);                                // result = result * 1.0/12.0;
    896    
    897    
     896
     897
    898898    if(n % 2 == 1) {
    899899        mpq_set_ui(qtemps, 1, 4);
    900900        mpq_sub(result, result, qtemps);                            // result = result - .25;
     
    910910        mpfr_set_ui(result, 1, round_mode);                         //result = 1
    911911        return;
    912912    }
    913    
     913
    914914    mpfr_set_ui(result, 0, round_mode);
    915915
    916916    unsigned int h = 0;
    917917    for(h = 1; h < k+1; h++) {
    918918        if(GCD(h,k) == 1) {
    919            
     919
    920920            // Note that we compute each term of the summand as
    921921            //      result += cos(pi * ( s(h,k) - (2.0 * h * n)/k) );
    922922            //
     
    924924            // down in the introduction, and we don't need to compute
    925925            // the imaginary part because we know that, in the end, the
    926926            // imaginary part will be 0, as we are computing an integer.
    927            
     927
    928928            q_s(qtempa, h, k);
    929              
     929
    930930            //mpfr_mul_q(tempa1, mp_pi, qtempa, round_mode);
    931931            //mpfr_mul_ui(tempa1, tempa1, k * k, round_mode);
    932932
    933933            //mpfr_set_q(tempa1, qtempa, round_mode);
    934             unsigned int r = n % k;                                     // here we make use of the fact that the 
     934            unsigned int r = n % k;                                     // here we make use of the fact that the
    935935            unsigned int d = GCD(r,k);                                  // cos() term written above only depends
    936936            unsigned int K;                                             // on {hn/k}.
    937937            if(d > 1) {
     
    949949            }
    950950            mpq_set_ui(qtempa2, h*r, K);
    951951            mpq_sub(qtempa, qtempa, qtempa2);
    952            
     952
    953953            //mpfr_set_q(tempa2, qtempa, round_mode);                   // This might be faster, according to
    954954            //cospi(tempa1, tempa2);                                    // the comments in Ralf Stephan's part.c, but
    955955                                                                        // I haven't noticed a significant speed up.
    956956                                                                        // (Perhaps a different version that takes an mpq_t
    957957                                                                        // as an input might be faster.)
    958            
     958
    959959            mpfr_mul_q(tempa1, mp_pi, qtempa, round_mode);
    960960            mpfr_cos(tempa1, tempa1, round_mode);
    961961            mpfr_add(result, result, tempa1, round_mode);
    962  
     962
    963963        }
    964        
     964
    965965    }
    966966
    967967}
    968968
    969 template <class T> 
     969template <class T>
    970970inline T partial_sum_of_t(unsigned int n, unsigned int &k, unsigned int exit_precision, unsigned int extra_precision, double error = 0) {
    971971    unsigned int current_precision = compute_current_precision(n, k - 1, extra_precision);
    972972    T result = 0;
     
    996996    // we can find p(n) by rounding.
    997997    //
    998998    //
    999     // NOTE: result should NOT have been initialized when this is called, 
     999    // NOTE: result should NOT have been initialized when this is called,
    10001000    // as we initialize it to the proper precision in this function.
    1001    
     1001
    10021002    double error = .25;
    10031003    int extra = compute_extra_precision(n, error);
    10041004
     
    10101010    mpfr_init2(t2, initial_precision);                              //
    10111011    mpfr_init2(result, initial_precision);                          //
    10121012    mpfr_set_ui(result, 0, round_mode);                             //
    1013                                                                    
     1013
    10141014    initialize_mpz_and_mpq_variables();
    10151015    initialize_constants(initial_precision, n);                     // Now that we have the precision information, we initialize some constants
    10161016                                                                    // that will be used throughout, and also
     
    10191019
    10201020    unsigned int current_precision = initial_precision;
    10211021    unsigned int new_precision;
    1022    
     1022
    10231023    // We start by computing with high precision arithmetic, until
    10241024    // we are sure enough that we don't need that much precision
    10251025    // anymore. Once current_precision == min_precision, we drop
     
    10301030
    10311031    unsigned int k = 1;                                             // (k holds the index of the summand that we are computing.)
    10321032    for(k = 1; current_precision > level_two_precision; k++) {      //
    1033        
     1033
    10341034        mpfr_sqrt_ui(t1, k, round_mode);                            // t1 = sqrt(k)
    10351035                                                                    //
    10361036        mp_a(t2, n, k);                                             // t2 = A_k(n)
    1037        
     1037
    10381038        if(debuga) {
    10391039            cout << "a(" << k <<  ") = ";
    10401040            mpfr_out_str(stdout, 10, 10, t2, round_mode);
    10411041            cout << endl;
    10421042        }
    1043        
     1043
    10441044        mpfr_mul(t1, t1, t2, round_mode);                           // t1 = sqrt(k)*A_k(n)
    10451045                                                                    //
    10461046        mp_f(t2, k);                                                // t2 = f_k(n)
     
    10891089    long double ld_partial_sum = partial_sum_of_t<long double>(n,k,level_five_precision, extra, 0);
    10901090    mpfr_set_ld(t1, ld_partial_sum, round_mode);
    10911091    mpfr_add(result, result, t1, round_mode);
    1092    
     1092
    10931093    double d_partial_sum = partial_sum_of_t<double>(n,k,0,extra,error);
    10941094
    10951095
     
    11041104    clear_mpz_and_mpq_variables();
    11051105    clear_mpfr_variables();
    11061106    mpfr_clear(t1);
    1107     mpfr_clear(t2);                       
     1107    mpfr_clear(t2);
    11081108}
    11091109
    11101110template <class T>
     
    11181118        return 1;
    11191119    }
    11201120    T result = 0;
    1121  
     1121
    11221122    unsigned int h = 0;
    11231123    for(h = 1; h < k+1; h++) {
    11241124        if(GCD(h,k) == 1) {
     
    11341134T s(unsigned int h, unsigned int k) {
    11351135    //mpq_set_ui(result, 1, 1);
    11361136    //return;
    1137    
     1137
    11381138    //return T(0);  // Uncommenting this line saves 25% or more time in the computation of p(10^9) in the current version (.51)
    11391139                    // but, of course, gives wrong results. (This is just to give an idea of how much time could
    11401140                    // possibly be saved by optimizing this function.
    1141    
     1141
    11421142    if(k < 3) {
    11431143        return T(0);
    11441144    }
     
    11861186
    11871187    //if(k % h == 2) {
    11881188    //}
    1189    
     1189
    11901190    int r1 = k;
    11911191    int r2 = h;
    11921192
     
    11991199                    // we only need to make sure that r2 > 0
    12001200    {
    12011201        temp = T(int(r1 * r1 + r2 * r2 + 1))/(n * r1 * r2);
    1202         temp3 = r1 % r2;                                           
    1203         r1 = r2;                                                   
    1204         r2 = temp3;                                                 
    1205        
     1202        temp3 = r1 % r2;
     1203        r1 = r2;
     1204        r2 = temp3;
     1205
    12061206        result += temp;                        // result += temp1;
    12071207
    12081208        n = -n;
    12091209    }
    12101210
    12111211    result *= one_over_12<T>();
    1212    
     1212
    12131213    if(n < 0) {
    12141214        result -= T(.25);
    12151215    }
     
    12371237      v = b;
    12381238      do {
    12391239         t = u % v;
    1240          u = v; 
     1240         u = v;
    12411241         v = t;
    12421242      } while (v != 0);
    12431243
     
    12521252// The following function was copied from Ralf Stephan's code,
    12531253// mentioned in the introduction, and then some variable names
    12541254// were changed. The function is not currently used, however.
    1255 void cospi (mpfr_t res, 
    1256                 mpfr_t x)
     1255void cospi (mpfr_t res,
     1256        mpfr_t x)
    12571257{
    1258 //      mpfr_t t, tt, half, fourth;
    1259        
    1260 //      mpfr_init2 (t, prec);
    1261 //      mpfr_init2 (tt, prec);
    1262 //      mpfr_init2 (half, prec);
    1263 //      mpfr_init2 (fourth, prec);
     1258//    mpfr_t t, tt, half, fourth;
    12641259
    1265 //      mpfr_set_ui (half, 1, r);
    1266 //      mpfr_div_2ui (half, half, 1, r);
    1267 //      mpfr_div_2ui (fourth, half, 1, r);
     1260//    mpfr_init2 (t, prec);
     1261//    mpfr_init2 (tt, prec);
     1262//    mpfr_init2 (half, prec);
     1263//    mpfr_init2 (fourth, prec);
     1264
     1265//    mpfr_set_ui (half, 1, r);
     1266//    mpfr_div_2ui (half, half, 1, r);
     1267//    mpfr_div_2ui (fourth, half, 1, r);
    12681268
    12691269        // NOTE: switched t to tempc2
    12701270        //       and tt to tempc1
    12711271
    12721272
    12731273        mp_rnd_t r = round_mode;
    1274        
    12751274
    1276         mpfr_div_2ui (tempc1, x, 1, r);
    1277         mpfr_floor (tempc1, tempc1);
    1278         mpfr_mul_2ui (tempc1, tempc1, 1, r);
    1279         mpfr_sub (tempc2, x, tempc1, r);
    1280         if (mpfr_cmp_ui (tempc2, 1) > 0)
    1281                 mpfr_sub_ui (tempc2, tempc2, 2, r);
    1282         mpfr_abs (tempc1, tempc2, r);
    1283         if (mpfr_cmp (tempc1, half) > 0)
    1284         {
    1285                 mpfr_ui_sub (tempc2, 1, tempc1, r);
    1286                 mpfr_abs (tempc1, tempc2, r);
    1287                 if (mpfr_cmp (tempc1, fourth) > 0)
    1288                 {
    1289                         if (mpfr_sgn (tempc2) > 0)
    1290                                 mpfr_sub (tempc2, half, tempc2, r);
    1291                         else
    1292                                 mpfr_add (tempc2, tempc2, half, r);
    1293                         mpfr_mul (tempc2, tempc2, mp_pi, r);
    1294                         mpfr_sin (tempc2, tempc2, r);
    1295                 }
    1296                 else
    1297                 {
    1298                         mpfr_mul (tempc2, tempc2, mp_pi, r);
    1299                         mpfr_cos (tempc2, tempc2, r);   
    1300                 }
    1301                 mpfr_neg (res, tempc2, r);
    1302         }
    1303         else
    1304         {
    1305                 mpfr_abs (tempc1, tempc2, r);
    1306                 if (mpfr_cmp (tempc1, fourth) > 0)
    1307                 {
    1308                         if (mpfr_sgn (tempc2) > 0)
    1309                                 mpfr_sub (tempc2, half, tempc2, r);
    1310                         else
    1311                                 mpfr_add (tempc2, tempc2, half, r);
    1312                         mpfr_mul (tempc2, tempc2, mp_pi, r);
    1313                         mpfr_sin (res, tempc2, r);
    1314                 }
    1315                 else
    1316                 {
    1317                         mpfr_mul (tempc2, tempc2, mp_pi, r);
    1318                         mpfr_cos (res, tempc2, r);
    1319                 }
    1320         }
    13211275
    1322 //      mpfr_clear (half);
    1323 //      mpfr_clear (fourth);
    1324 //      mpfr_clear (t);
    1325 //      mpfr_clear (tt);
     1276    mpfr_div_2ui (tempc1, x, 1, r);
     1277    mpfr_floor (tempc1, tempc1);
     1278    mpfr_mul_2ui (tempc1, tempc1, 1, r);
     1279    mpfr_sub (tempc2, x, tempc1, r);
     1280    if (mpfr_cmp_ui (tempc2, 1) > 0)
     1281        mpfr_sub_ui (tempc2, tempc2, 2, r);
     1282    mpfr_abs (tempc1, tempc2, r);
     1283    if (mpfr_cmp (tempc1, half) > 0)
     1284    {
     1285        mpfr_ui_sub (tempc2, 1, tempc1, r);
     1286        mpfr_abs (tempc1, tempc2, r);
     1287        if (mpfr_cmp (tempc1, fourth) > 0)
     1288        {
     1289            if (mpfr_sgn (tempc2) > 0)
     1290                mpfr_sub (tempc2, half, tempc2, r);
     1291            else
     1292                mpfr_add (tempc2, tempc2, half, r);
     1293            mpfr_mul (tempc2, tempc2, mp_pi, r);
     1294            mpfr_sin (tempc2, tempc2, r);
     1295        }
     1296        else
     1297        {
     1298            mpfr_mul (tempc2, tempc2, mp_pi, r);
     1299            mpfr_cos (tempc2, tempc2, r);
     1300        }
     1301        mpfr_neg (res, tempc2, r);
     1302    }
     1303    else
     1304    {
     1305        mpfr_abs (tempc1, tempc2, r);
     1306        if (mpfr_cmp (tempc1, fourth) > 0)
     1307        {
     1308            if (mpfr_sgn (tempc2) > 0)
     1309                mpfr_sub (tempc2, half, tempc2, r);
     1310            else
     1311                mpfr_add (tempc2, tempc2, half, r);
     1312            mpfr_mul (tempc2, tempc2, mp_pi, r);
     1313            mpfr_sin (res, tempc2, r);
     1314        }
     1315        else
     1316        {
     1317            mpfr_mul (tempc2, tempc2, mp_pi, r);
     1318            mpfr_cos (res, tempc2, r);
     1319        }
     1320    }
     1321
     1322//    mpfr_clear (half);
     1323//    mpfr_clear (fourth);
     1324//    mpfr_clear (t);
     1325//    mpfr_clear (tt);
    13261326}
    13271327
    13281328/* answer must have already been mpz_init'd. */
     
    13321332        return 0;
    13331333    }
    13341334    mpfr_t result;
    1335    
     1335
    13361336    mp_t(result, n);
    1337    
     1337
    13381338    mpfr_get_z(answer, result, round_mode);
    13391339
    13401340    mpfr_clear(result);
     
    13441344
    13451345int main(int argc, char *argv[]){
    13461346    //init();
    1347    
     1347
    13481348    unsigned int n = 10;
    13491349
    13501350    if(argc > 1)
     
    13831383        return 0;
    13841384    }
    13851385    //mpfr_t result;
    1386    
     1386
    13871387    //mp_t(result, n);
    1388    
     1388
    13891389    mpz_t answer;
    13901390    mpz_init(answer);
    13911391    part(answer, n);
     
    13931393    //mpfr_get_z(answer, result, round_mode);
    13941394
    13951395    mpz_out_str (stdout, 10, answer);
    1396    
     1396
    13971397    cout << endl;
    13981398
    13991399    return 0;
     
    14121412
    14131413    mpz_t expected_value;
    14141414    mpz_t actual_value;
    1415    
     1415
    14161416    mpz_init(expected_value);
    14171417    mpz_init(actual_value);
    14181418
     
    14351435
    14361436    if(mpz_cmp(expected_value, actual_value) != 0)
    14371437        return n;
    1438    
     1438
    14391439    cout << " OK." << endl;
    14401440
    14411441    n = 100;
     
    14791479
    14801480    if(mpz_cmp(expected_value, actual_value) != 0)
    14811481        return n;
    1482    
     1482
    14831483    cout << " OK." << endl;
    14841484
    14851485    n = 1000000;
     
    14901490
    14911491    if(mpz_cmp(expected_value, actual_value) != 0)
    14921492        return n;
    1493    
     1493
    14941494    cout << " OK." << endl;
    14951495
    14961496
    14971497    // We now run some tests based on the fact that if n = 369 (mod 385) then p(n) = 0 (mod 385).
    1498    
     1498
    14991499    n = 369 + 10*385;
    15001500    cout << "Computing p(" << n << ")...";
    15011501    cout.flush();
     
    15531553    // Randomized testing
    15541554
    15551555    srand( time(NULL) );
    1556    
     1556
    15571557    for(int i = 0; i < 100; i++) {
    15581558        n = int(100000 * double(rand())/double(RAND_MAX) + 1);
    15591559        n = n - (n % 385) + 369;
     
    16491649        part(actual_value, n);
    16501650        if(mpz_divisible_ui_p(actual_value, 385) == 0)
    16511651            return n;
    1652        
     1652
    16531653        cout << " OK." << endl;
    16541654
    16551655        for(int i = 0; i < 10; i++) {
     
    16631663            }
    16641664            cout << " OK." << endl;
    16651665        }
    1666        
     1666
    16671667        n = 1000000000 + 139;
    16681668        cout << "Computing p(" << n << ")...";
    16691669        cout.flush();
    16701670        part(actual_value, n);
    16711671        if(mpz_divisible_ui_p(actual_value, 385) == 0)
    16721672            return n;
    1673        
     1673
    16741674        cout << " OK." << endl;
    16751675
    16761676
  • sage/geometry/triangulation/data.cc

    diff --git a/sage/geometry/triangulation/data.cc b/sage/geometry/triangulation/data.cc
    a b  
    1212
    1313// ------- class vertices ----------------------------
    1414
    15 int vertices::n=0; 
     15int vertices::n=0;
    1616
    17 int vertices::d=0; 
     17int vertices::d=0;
    1818
    1919vertices_lookup vertices::lookup=vertices_lookup();
    2020
    2121
    22 void vertices::set_dimensions(int N, int D) 
    23 { 
    24   n=N; d=D; 
     22void vertices::set_dimensions(int N, int D)
     23{
     24  n=N; d=D;
    2525  lookup.generate_tables(n,d);
    2626}
    2727
     
    4141{
    4242  simplex s=1;
    4343  vertex k=1, l;
    44   const_iterator l_iterator=begin(); 
    45   for (vertex i=1; i<=d; ++i) {   
    46     l = (*l_iterator)+1;   
     44  const_iterator l_iterator=begin();
     45  for (vertex i=1; i<=d; ++i) {
     46    l = (*l_iterator)+1;
    4747    for (vertex j=k; j<=l-1; ++j)
    48       s=s+lookup.get_binomial(n-j,d-i);   
    49     k=l+1;   
    50     ++l_iterator; 
    51   } 
     48      s=s+lookup.get_binomial(n-j,d-i);
     49    k=l+1;
     50    ++l_iterator;
     51  }
    5252  return(s);
    5353}
    5454
    5555
    56 const vertices & vertices::simplex_to_vertices(const simplex & S) 
     56const vertices & vertices::simplex_to_vertices(const simplex & S)
    5757{
    5858  (*this) = lookup.simplex_to_vertices(S);
    59   return(*this); 
     59  return(*this);
    6060}
    6161
    6262
    6363std::ostream & operator << (std::ostream & out, const vertices & v)
    6464{
    65   out << *( v.begin() ); 
    66   vertices::const_iterator i=v.begin();   i++; 
     65  out << *( v.begin() );
     66  vertices::const_iterator i=v.begin();   i++;
    6767  for (; i!=v.end(); ++i)
    6868    out << "_" << *i;
    6969  return(out);
     
    7777  if (n1>n2) return(false);
    7878  // a,b have same size
    7979  for (vertices::iterator i=a.begin(), j=b.begin();
    80        i!=a.end() && j!=b.end(); ++i, ++j) {   
    81     if (*i < *j) return(true); 
     80       i!=a.end() && j!=b.end(); ++i, ++j) {
     81    if (*i < *j) return(true);
    8282    if (*i > *j) return(false);
    8383  }
    8484  // a == b
     
    8686}
    8787
    8888
    89 bool operator==(const vertices & a, const vertices & b) 
     89bool operator==(const vertices & a, const vertices & b)
    9090{
    91   return( std::set<vertex,std::less<vertex> >(a) == 
     91  return( std::set<vertex,std::less<vertex> >(a) ==
    9292          std::set<vertex,std::less<vertex> >(b) );
    9393}
    9494
    9595
    9696// ------- class vertices_lookup --------------------
    9797
    98 const vertices & vertices_lookup::simplex_to_vertices(const simplex & S) const 
     98const vertices & vertices_lookup::simplex_to_vertices(const simplex & S) const
    9999{
    100100  return( SimplexToVertices[S-1] );
    101101}
     
    104104// this one is messy, but it only reverses vertices_to_simplex()
    105105vertices vertices_lookup::manual_vertices_to_simplex(const simplex & S) const
    106106{
    107   vertices result; 
     107  vertices result;
    108108  simplex s=S;
    109   simplex b; 
     109  simplex b;
    110110  vertex i,j,l=0,k;
    111   for (k=1; k<d; k++) {   
    112     l++;  i=l;   j=1; 
     111  for (k=1; k<d; k++) {
     112    l++;  i=l;   j=1;
    113113    b=binomial(n-l,d-k);
    114     while (s>b && b>0) { 
     114    while (s>b && b>0) {
    115115      j++;  l++;
    116       s=s-b; 
     116      s=s-b;
    117117      b=binomial(n-l,d-k);
    118118    }
    119119    result.insert(result.begin(),l-1);
    120120  }
    121121  result.insert(result.begin(),s+l-1);
    122122  return(result);
    123 } 
     123}
    124124
    125 // tweaked for special (e.g. negative) values 
    126 int vertices_lookup::get_binomial(int i, int j) const 
    127 { 
     125// tweaked for special (e.g. negative) values
     126int vertices_lookup::get_binomial(int i, int j) const
     127{
    128128  if (i>=0 && i<=n && j>=0 && j<=d && j<=i)
    129     return(fast_binomial[i][j]); 
    130   else 
     129    return(fast_binomial[i][j]);
     130  else
    131131    return(1);
    132 } 
     132}
    133133
    134134
    135135void vertices_lookup::generate_tables(int N, int D)
    136136{
    137   n=N;  d=D; 
     137  n=N;  d=D;
    138138  fast_binomial.erase(fast_binomial.begin(),fast_binomial.end());
    139139  for(int i=0; i<=n; ++i) {
    140     std::vector<int> row;   
     140    std::vector<int> row;
    141141    for(int j=0; j<=i && j<=d; j++)
    142142      row.push_back( binomial(i,j) );
    143143    fast_binomial.push_back( row );
    144144  }
    145   SimplexToVertices.erase(SimplexToVertices.begin(),SimplexToVertices.end()); 
     145  SimplexToVertices.erase(SimplexToVertices.begin(),SimplexToVertices.end());
    146146  for(int s=1; s<=binomial(n,d); ++s)
    147147    SimplexToVertices.push_back(manual_vertices_to_simplex(s));
    148148}
     
    152152
    153153// ------- class compact_simplices -------------------
    154154
    155 compact_simplices::compact_simplices() 
    156 { 
    157   reserve(50); 
     155compact_simplices::compact_simplices()
     156{
     157  reserve(50);
    158158}
    159  
     159
    160160compact_simplices::~compact_simplices()
    161161{
    162162}
    163  
     163
    164164std::ostream & operator << (std::ostream & out, const compact_simplices & t)
    165165{
    166   out << "[" << *( t.begin() ); 
     166  out << "[" << *( t.begin() );
    167167  for (compact_simplices::const_iterator i=t.begin()+1; i!=t.end(); ++i)
    168168    out << "," << *i;
    169169  out << "]";
     
    179179hash_value compact_simplices::hash_function() const
    180180{
    181181  hash_value result=0;
    182   int n=101; 
     182  int n=101;
    183183  for (const_iterator i=begin(); i!=end(); ++i, n+=37)
    184184    result += n* (*i) + n*n;
    185185  return(result);
     
    201201{
    202202  v.erase(v.begin(), v.end());
    203203  copy(s.begin(), s.end(), back_inserter(v) );
    204   compress(); 
     204  compress();
    205205}
    206206
    207207
     
    211211
    212212void simplices::compress()
    213213{
    214   erase(begin(),end()); 
     214  erase(begin(),end());
    215215  for (std::vector<vertices>::const_iterator i=v.begin(); i!=v.end(); i++)
    216     push_back( (*i).vertices_to_simplex() ); 
    217   std::sort(begin(),end()); 
     216    push_back( (*i).vertices_to_simplex() );
     217  std::sort(begin(),end());
    218218}
    219219
    220220
    221221void simplices::decompress()
    222222{
    223   v.erase(v.begin(),v.end()); 
     223  v.erase(v.begin(),v.end());
    224224  for (const_iterator i=begin(); i!=end(); i++)
    225225    v.push_back( vertices().simplex_to_vertices(*i) );
    226226}
     
    228228
    229229bool simplices::starshaped(const vertex origin) const
    230230{
    231   bool result = true; 
     231  bool result = true;
    232232  for (std::vector<vertices>::const_iterator
    233         vi=v.begin(); vi!=v.end(); vi++) {
     233        vi=v.begin(); vi!=v.end(); vi++) {
    234234    result=result && (find(vi->begin(),vi->end(),origin)!=vi->end());
    235235  }
    236   return(result);           
     236  return(result);
    237237}
    238238
    239239
     
    241241{
    242242  vertices support;
    243243  for (std::vector<vertices>::const_iterator
    244         vi=v.begin(); vi!=v.end(); vi++) {
     244        vi=v.begin(); vi!=v.end(); vi++) {
    245245    support.insert(vi->begin(), vi->end());
    246246  }
    247   return support.full_set(); 
     247  return support.full_set();
    248248}
    249249
    250250
    251251std::ostream & operator << (std::ostream & out, const simplices & s)
    252252{
    253   out << "[" << s.v.front(); 
    254   for (std::vector<vertices>::const_iterator 
    255         si=s.v.begin()+1; si!=s.v.end(); si++)
     253  out << "[" << s.v.front();
     254  for (std::vector<vertices>::const_iterator
     255        si=s.v.begin()+1; si!=s.v.end(); si++)
    256256    out << ", " << *si;
    257257  out << "]";
    258258  return(out);
     
    265265  deltaplus.reserve(10);
    266266  deltaminus.reserve(10);
    267267}
    268  
     268
    269269flip::flip(const std::vector<vertices>& pos, const std::vector<vertices>& neg)
    270270  : deltaplus(pos), deltaminus(neg)
    271271{}
     
    274274{
    275275  deltaplus.reserve(10);
    276276  deltaminus.reserve(10);
    277   (*this)=copyfrom; 
     277  (*this)=copyfrom;
    278278}
    279279
    280280flip::~flip()
     
    284284flip & flip::operator =(const flip & copyfrom)
    285285{
    286286  deltaplus  = copyfrom.get_deltaplus();
    287   deltaminus = copyfrom.get_deltaminus(); 
    288   return(*this); 
     287  deltaminus = copyfrom.get_deltaminus();
     288  return(*this);
    289289}
    290290
    291291std::ostream & operator << (std::ostream & out, const flip & f)
    292292{
    293   out << "< "; 
    294   for (std::vector<vertices>::const_iterator i=f.get_deltaplus().begin(); 
    295        i!=f.get_deltaplus().end(); ++i) {     
     293  out << "< ";
     294  for (std::vector<vertices>::const_iterator i=f.get_deltaplus().begin();
     295       i!=f.get_deltaplus().end(); ++i) {
    296296    std::cout << *i << " ";
    297   }     
    298   out << "| "; 
    299   for (std::vector<vertices>::const_iterator i=f.get_deltaminus().begin(); 
    300        i!=f.get_deltaminus().end(); ++i) { 
    301     std::cout << *i << " ";   
    302   }     
     297  }
     298  out << "| ";
     299  for (std::vector<vertices>::const_iterator i=f.get_deltaminus().begin();
     300       i!=f.get_deltaminus().end(); ++i) {
     301    std::cout << *i << " ";
     302  }
    303303  out << ">";
    304304  return(out);
    305305}
    306306
    307307void flip::mirror()
    308308{
    309   swap(deltaplus,deltaminus); 
     309  swap(deltaplus,deltaminus);
    310310}
    311311
    312312
     
    328328
    329329goodcircuit::goodcircuit(const simplices & s, const flip & f)
    330330  :supporter(f)
    331 {   
    332   const std::vector<vertices> & deltaplus = f.get_deltaplus(); 
    333   const std::vector<vertices> & v = s.get_vertices(); 
     331{
     332  const std::vector<vertices> & deltaplus = f.get_deltaplus();
     333  const std::vector<vertices> & v = s.get_vertices();
    334334  supported.reserve(deltaplus.size());
    335   good=true; 
    336  
     335  good=true;
    337336
    338   // find simplices that include members of deltaplus
    339   std::vector<vertices> empty_result; 
     337
     338  // find simplices that include members of deltaplus
     339  std::vector<vertices> empty_result;
    340340  for (size_t n=0; n<deltaplus.size(); n++) {
    341     supported.push_back(empty_result);  supported[n].reserve(10); 
    342     bool found_match=false;   
     341    supported.push_back(empty_result);  supported[n].reserve(10);
     342    bool found_match=false;
    343343    for (std::vector<vertices>::const_iterator i=v.begin(); i!=v.end(); ++i) {
    344       if (includes( (*i).begin(), (*i).end(), 
    345                     deltaplus[n].begin(), deltaplus[n].end() )) {
    346         found_match=true;
    347         supported[n].push_back(*i);     
     344      if (includes( (*i).begin(), (*i).end(),
     345                    deltaplus[n].begin(), deltaplus[n].end() )) {
     346        found_match=true;
     347        supported[n].push_back(*i);
    348348      }
    349349    }
    350350    good=good && found_match;
    351     if (!good) return;   
     351    if (!good) return;
    352352  }
    353353
    354354  // check that all links are equal
    355355  for (size_t n=0; n<deltaplus.size(); n++) {
    356     std::set<vertices,vertices_order> l; 
     356    std::set<vertices,vertices_order> l;
    357357    for (std::vector<vertices>::const_iterator i=supported[n].begin();
    358          i!=supported[n].end(); ++i) {     
    359       vertices one_simplex;   
    360       std::insert_iterator<vertices> 
    361         ins(one_simplex,one_simplex.begin());   
    362       set_difference( (*i).begin(), (*i).end(), 
    363                       deltaplus[n].begin(), deltaplus[n].end(),
    364                       ins );
    365       l.insert(l.begin(),one_simplex);     
    366     }   
    367     link.push_back(l);   
     358             i!=supported[n].end(); ++i) {
     359      vertices one_simplex;
     360      std::insert_iterator<vertices>
     361        ins(one_simplex,one_simplex.begin());
     362      set_difference( (*i).begin(), (*i).end(),
     363                       deltaplus[n].begin(), deltaplus[n].end(),
     364                       ins );
     365      l.insert(l.begin(),one_simplex);
     366    }
     367    link.push_back(l);
    368368    if (n>0 && link[0]!=link[n]) {
    369369      good=false;
    370370      return;
    371     }                                       
     371    }
    372372  }
    373373}
    374374
     
    376376void goodcircuit::do_flip(const simplices & s, const flip & f)
    377377{
    378378  bistellarneighbor.erase(bistellarneighbor.begin(), bistellarneighbor.end());
    379  
     379
    380380  const std::set<vertices,vertices_order> & use_link = link[0];
    381   const std::vector<vertices> & v = s.get_vertices(); 
    382   const std::vector<vertices> & deltaplus  = f.get_deltaplus(); 
    383   const std::vector<vertices> & deltaminus = f.get_deltaminus(); 
    384  
    385   // start with all simplices 
    386   std::set<vertices,vertices_order> tri; 
     381  const std::vector<vertices> & v = s.get_vertices();
     382  const std::vector<vertices> & deltaplus  = f.get_deltaplus();
     383  const std::vector<vertices> & deltaminus = f.get_deltaminus();
     384
     385  // start with all simplices
     386  std::set<vertices,vertices_order> tri;
    387387  copy(v.begin(), v.end(), inserter(tri,tri.begin()) );
    388  
     388
    389389  // remove all suported simplices
    390390  std::set<vertices,vertices_order> to_remove;
    391  
     391
    392392  for (std::set<vertices,vertices_order>::const_iterator i=use_link.begin();
    393        i!=use_link.end(); ++i) 
     393       i!=use_link.end(); ++i)
    394394    for (std::vector<vertices>::const_iterator j=deltaplus.begin();
    395          j!=deltaplus.end(); ++j) {   
     395         j!=deltaplus.end(); ++j) {
    396396      vertices one_simplex;
    397397      set_union( (*j).begin(), (*j).end(),
    398                 (*i).begin(), (*i).end(),
    399                 inserter(one_simplex,one_simplex.begin()) );
    400       to_remove.insert(to_remove.begin(),one_simplex);           
     398                (*i).begin(), (*i).end(),
     399                inserter(one_simplex,one_simplex.begin()) );
     400      to_remove.insert(to_remove.begin(),one_simplex);
    401401    }
    402   set_difference( tri.begin(), tri.end(), 
    403                   to_remove.begin(), to_remove.end(),
    404                   inserter(bistellarneighbor,bistellarneighbor.begin()),
    405                   vertices_order() );
    406    
     402  set_difference( tri.begin(), tri.end(),
     403                  to_remove.begin(), to_remove.end(),
     404                  inserter(bistellarneighbor,bistellarneighbor.begin()),
     405                  vertices_order() );
     406
    407407  // now add the flip of the supported simplices
    408408  for (std::set<vertices,vertices_order>::const_iterator i=use_link.begin();
    409409       i!=use_link.end(); ++i)
    410410    for (std::vector<vertices>::const_iterator j=deltaminus.begin();
    411         j!=deltaminus.end(); j++) {
     411        j!=deltaminus.end(); j++) {
    412412      vertices one_simplex;
    413413      set_union( (*j).begin(), (*j).end(),
    414                 (*i).begin(), (*i).end(),
    415                 inserter(one_simplex,one_simplex.begin()) );
     414                (*i).begin(), (*i).end(),
     415                inserter(one_simplex,one_simplex.begin()) );
    416416      bistellarneighbor.insert( bistellarneighbor.begin(), one_simplex );
    417     } 
     417    }
    418418}
    419419
    420420
  • sage/geometry/triangulation/triangulations.cc

    diff --git a/sage/geometry/triangulation/triangulations.cc b/sage/geometry/triangulation/triangulations.cc
    a b  
    1111    star(-1),
    1212    fine(false),
    1313    need_resize(false)
    14 {} 
     14{}
    1515
    1616
    1717// find the correct position for t in the hash
    18 void triangulations::find_hash_position(const compact_simplices& t, 
    19                                         hash_value& pos, bool& is_new) const
     18void triangulations::find_hash_position(const compact_simplices& t,
     19                                        hash_value& pos, bool& is_new) const
    2020{
    2121  hash_value freespace;
    2222  const hash_value initial_guess = t.hash_function() % hash_max;
    2323
    2424  for (hash_value i=0; i<hash_max; ++i) {
    2525    pos = (initial_guess+i) % hash_max;
    26     if (hash_list[pos]==hash_max) { 
     26    if (hash_list[pos]==hash_max) {
    2727      // found empty place in hash
    2828      is_new=true;
    2929      if (i>5)
    30         need_resize = true;
     30        need_resize = true;
    3131      return;
    32     } 
     32    }
    3333    else if ((*this)[hash_list[pos]]==t) {
    3434      // found ourselves in hash
    3535      is_new = false;
     
    5454    for (size_t i=0; i<size(); i++) {
    5555      find_hash_position( (*this)[i], pos, is_new);
    5656      assert(is_new);
    57       hash_list[pos] = i; 
     57      hash_list[pos] = i;
    5858    }
    5959    need_resize = false;
    6060    find_hash_position(new_triang, pos, is_new);
    6161  }
    6262
    63   push_back(new_triang);   
    64   hash_list[pos] = size()-1; 
     63  push_back(new_triang);
     64  hash_list[pos] = size()-1;
    6565}
    6666
    6767
    6868void triangulations::add_neighbours(const simplices & s)
    6969{
    70   for (flips::const_iterator 
    71         f=bistellar_flips.begin(); f!=bistellar_flips.end(); ++f) {
    72     goodcircuit goody(s,*f);   
     70  for (flips::const_iterator
     71        f=bistellar_flips.begin(); f!=bistellar_flips.end(); ++f) {
     72    goodcircuit goody(s,*f);
    7373    if (goody.is_good()) {
    7474      goody.do_flip(s,*f);
    75       compact_simplices new_triang=goody.get_neighbor();   
    76       add_triang_if_new(new_triang);     
    77     }   
    78   } 
     75      compact_simplices new_triang=goody.get_neighbor();
     76      add_triang_if_new(new_triang);
     77    }
     78  }
    7979}
    8080
    8181
     
    8484  while (position != this->size()) {
    8585     // eat all non-star triangulations
    8686    simplices triangulation((*this)[position]);
    87    
     87
    8888    if ((star>=0) && !(triangulation.starshaped(star))) {
    8989      next_triangulation();
    9090      continue;
     
    113113triangulations_ptr init_triangulations
    114114(int n, int d, int star, bool fine, PyObject* py_seed, PyObject* py_flips)
    115115{
    116   vertices().set_dimensions(n,d); 
     116  vertices().set_dimensions(n,d);
    117117
    118118  compact_simplices seed;
    119119  for (int i=0; i<PySequence_Size(py_seed); i++) {
     
    133133      PyObject* py_simplex = PySequence_GetItem(py_flip_pos,j);
    134134      vertices simplex;
    135135      for (int k=0; k<PySequence_Size(py_simplex); k++) {
    136         PyObject* py_vertex = PySequence_GetItem(py_simplex,k);
    137         simplex.insert(simplex.begin(), PyInt_AS_LONG(py_vertex));
    138         Py_DECREF(py_vertex);
     136        PyObject* py_vertex = PySequence_GetItem(py_simplex,k);
     137        simplex.insert(simplex.begin(), PyInt_AS_LONG(py_vertex));
     138        Py_DECREF(py_vertex);
    139139      }
    140140      pos.push_back(simplex);
    141141      Py_DECREF(py_simplex);
     
    146146      PyObject* py_simplex = PySequence_GetItem(py_flip_neg,j);
    147147      vertices simplex;
    148148      for (int k=0; k<PySequence_Size(py_simplex); k++) {
    149         PyObject* py_vertex = PySequence_GetItem(py_simplex,k);
    150         simplex.insert(simplex.begin(), PyInt_AS_LONG(py_vertex));
    151         Py_DECREF(py_vertex);
     149        PyObject* py_vertex = PySequence_GetItem(py_simplex,k);
     150        simplex.insert(simplex.begin(), PyInt_AS_LONG(py_vertex));
     151        Py_DECREF(py_vertex);
    152152      }
    153153      neg.push_back(simplex);
    154154      Py_DECREF(py_simplex);
     
    164164  // print data so we can see that it worked
    165165  /*
    166166  std::cout << "\n" << "Seed triangulation: " << seed << "\n";
    167   {   
    168     std::cout << "Vertices of seed triangulation: ";   
    169     simplices s(seed);   
    170     std::cout << s << "\n";   
    171   }       
    172   std::cout << "All flips:\n"; 
     167  {
     168    std::cout << "Vertices of seed triangulation: ";
     169    simplices s(seed);
     170    std::cout << s << "\n";
     171  }
     172  std::cout << "All flips:\n";
    173173  for (flips::const_iterator i=all_flips.begin(); i!=all_flips.end(); ++i)
    174174    std::cout << *i << std::endl;
    175175  */
     
    182182  if (fine)
    183183    t->require_fine_triangulation(fine);
    184184
    185   t->add_triang_if_new(seed); 
     185  t->add_triang_if_new(seed);
    186186
    187187  return t;
    188188}
     
    197197  PyObject* py_triang = PyTuple_New(triang.size());
    198198  for (int i=0; i<triang.size(); i++)
    199199    PyTuple_SET_ITEM(py_triang, i, PyInt_FromLong(triang[i]));
    200  
     200
    201201  return py_triang;
    202202}
    203203
  • sage/graphs/modular_decomposition/src/dm.c

    diff --git a/sage/graphs/modular_decomposition/src/dm.c b/sage/graphs/modular_decomposition/src/dm.c
    a b  
    3838#include <stdio.h>
    3939#include <stdlib.h>
    4040
    41 #define DEBUG 0                 /* si 0 aucune sortie graphique!! */
     41#define DEBUG 0                        /* si 0 aucune sortie graphique!! */
    4242
    43 /* dm.h definit les constantes FEUILLE, MODULE, etc... 
     43/* dm.h definit les constantes FEUILLE, MODULE, etc...
    4444ainsi que les structures noeud et fils. Les autres
    4545structures n'ont pas a etre connues par les programmes
    4646exterieurs et sont donc definies ici. */
    4747
    4848
    49 /* un sommet du graphe 
     49/* un sommet du graphe
    5050(utilise dans la premiere partie seulement, ainsi que ce qui suit)*/
    5151typedef struct Sommet {
    5252  int place;/* numero du sommet dans la NOUVELLE numerotation */
     
    6969    int debut;
    7070    int fin;
    7171    struct Sommet *firstpivot;
    72     int inpivot;                /*indice de la classe dans le tableau pivot */
    73     int inmodule;               /* (resp module); -1 si non present */
    74     int whereXa;                /* lie le couple X/Xa: vaut
    75                                    0 si X n'est actuellement lie a aucun Xa
    76                                    -1 si Xa est a gauche
    77                                    +1 si Xa est a droite */
    78     struct SClasse *suiv;       /* forment une liste chainee */
    79     struct SClasse *prec;       /*...doublement */
     72    int inpivot;                /*indice de la classe dans le tableau pivot */
     73    int inmodule;                /* (resp module); -1 si non present */
     74    int whereXa;                /* lie le couple X/Xa: vaut
     75                                   0 si X n'est actuellement lie a aucun Xa
     76                                   -1 si Xa est a gauche
     77                                   +1 si Xa est a droite */
     78    struct SClasse *suiv;        /* forment une liste chainee */
     79    struct SClasse *prec;        /*...doublement */
    8080} sclasse;
    8181
    8282/* plein de parametres statiques que algo1() donne a Raffine() */
     
    8989  int *n;
    9090} info;
    9191
    92 /* clef a deux entrees utilisee pour le tri lineaire 
     92/* clef a deux entrees utilisee pour le tri lineaire
    9393   represente l'arrete ij */
    9494typedef struct Clef2{
    9595  int i; //sommet pointeur
     
    123123    return (a > b) ? a : b;
    124124}
    125125/**************************************************************
    126 Premiere passe de l'algorithme: il s'agit de trouver une 
     126Premiere passe de l'algorithme: il s'agit de trouver une
    127127permutation factorisante du graphe. Nous utilisons les
    128128techniques de raffinement de partition. Tout cela est
    129129explique dans l'article de Habib, Viennot & Paul, dont je ne
     
    158158  nouv->inmodule = -1;
    159159  nouv->firstpivot = NULL;
    160160  nouv->prec = un;
    161   if (un != NULL)               /* accroche pas en tete de chaine */
     161  if (un != NULL)                /* accroche pas en tete de chaine */
    162162    un->suiv = nouv;
    163163  nouv->suiv = deux;
    164   if (deux != NULL)             /* pas en queue de chaine */
     164  if (deux != NULL)                /* pas en queue de chaine */
    165165    deux->prec = nouv;
    166166
    167167    /* debut et fin ne sont PAS initialises! */
     
    183183void Raffiner(sommet ** S, sommet * p, sommet * centre, info * I)
    184184{
    185185  /* melange raffiner, pivotset, insertright et addpivot */
    186   sadj *a;                      /* parcours l'adjacence du pivot */
    187   sommet *x;                    /* sommet quiva changer de classe */
    188   sclasse *X, *Xa;              /* x in X; Xa nouv classe de x */
     186  sadj *a;                        /* parcours l'adjacence du pivot */
     187  sommet *x;                        /* sommet quiva changer de classe */
     188  sclasse *X, *Xa;                /* x in X; Xa nouv classe de x */
    189189  sclasse *Z;
    190190  sclasse **pivot;
    191191  sclasse **module;
     
    204204    x = a->pointe;
    205205    X = x->classe;
    206206    if (X == p->classe)
    207       continue;         /* on raffine pas la classe du pivot! */
     207      continue;                /* on raffine pas la classe du pivot! */
    208208
    209209    if (X->whereXa == 0) {
    210210      /* c'est la premiere fois qu'on trouve un x
    211         appartenant a X lors de cet appel a raffiner */
     211        appartenant a X lors de cet appel a raffiner */
    212212
    213213      if ((centre->place < x->place && x->place < p->place)
    214           || (p->place < x->place && x->place < centre->place)) {
    215         /* insere a gauche */
    216         Xa = nouvclasse(X->prec, X);
    217         (*numclasse)++;
    218         permute(S, x->place, X->debut);
    219         X->debut++;
    220         X->whereXa = -1;
    221         Xa->whereXa = 1;        /* besoin dans le second tour */
     214          || (p->place < x->place && x->place < centre->place)) {
     215        /* insere a gauche */
     216        Xa = nouvclasse(X->prec, X);
     217        (*numclasse)++;
     218        permute(S, x->place, X->debut);
     219        X->debut++;
     220        X->whereXa = -1;
     221        Xa->whereXa = 1;        /* besoin dans le second tour */
    222222      }
    223       else {            /* insere a droite */
     223      else {                /* insere a droite */
    224224
    225         Xa = nouvclasse(X, X->suiv);
    226         (*numclasse)++;
    227         permute(S, x->place, X->fin);
    228         X->fin--;
    229         X->whereXa = 1;
    230         Xa->whereXa = -1;
     225        Xa = nouvclasse(X, X->suiv);
     226        (*numclasse)++;
     227        permute(S, x->place, X->fin);
     228        X->fin--;
     229        X->whereXa = 1;
     230        Xa->whereXa = -1;
    231231      }
    232232      x->classe = Xa;
    233233      Xa->debut = x->place;
     
    235235    }
    236236    else {
    237237      if (X->whereXa == -1) {
    238         Xa = X->prec;
    239         permute(S, x->place, X->debut);
    240         X->debut++;
    241         Xa->fin++;
     238        Xa = X->prec;
     239        permute(S, x->place, X->debut);
     240        X->debut++;
     241        Xa->fin++;
    242242      }
    243243      else {
    244         Xa = X->suiv;
    245         permute(S, x->place, X->fin);
    246         X->fin--;
    247         Xa->debut--;
     244        Xa = X->suiv;
     245        permute(S, x->place, X->fin);
     246        X->fin--;
     247        Xa->debut--;
    248248      }
    249249      x->classe = Xa;
    250250    }
     
    252252
    253253  for (a = p->adj; a != NULL; a = a->suiv)
    254254    /* deuxieme couche! Maintenant on va faire les addpivot,
    255            et remettre les whereXa a 0
    256            Noter qu'on lit les Xa et plus les X */
     255           et remettre les whereXa a 0
     256           Noter qu'on lit les Xa et plus les X */
    257257    {
    258258      x = a->pointe;
    259259      Xa = x->classe;
    260260      if (Xa->whereXa == 0)
    261         continue;               /* deja remis a zero! */
     261        continue;                /* deja remis a zero! */
    262262      if (Xa->whereXa == -1)
    263         X = Xa->prec;
     263        X = Xa->prec;
    264264      else
    265         X = Xa->suiv;
     265        X = Xa->suiv;
    266266
    267267      if (X->debut > X->fin) {
    268         /*on a trop enleve! X est vide
    269           -> on va le supprimer mechamment */
    270        
    271         (*numclasse)--;
    272         if (Xa->whereXa == 1) { /*deconnecte */
    273           Xa->suiv = X->suiv;
    274           if (Xa->suiv != NULL)
    275             Xa->suiv->prec = Xa;
    276         }
    277         else {
    278           Xa->prec = X->prec;
    279           if (Xa->prec != NULL)
    280             Xa->prec->suiv = Xa;
    281         }
    282         Xa->inpivot = X->inpivot;
    283         if (X->inpivot != -1)   /* ecrase X dans pivot */
    284           pivot[X->inpivot] = Xa;
    285         Xa->inmodule = X->inmodule;
    286         if (X->inmodule != -1)  /* ecrase X dans pivot */
    287           module[X->inmodule] = Xa;
     268        /*on a trop enleve! X est vide
     269          -> on va le supprimer mechamment */
    288270
    289         Xa->whereXa = 0;
    290         continue;
     271        (*numclasse)--;
     272        if (Xa->whereXa == 1) {        /*deconnecte */
     273          Xa->suiv = X->suiv;
     274          if (Xa->suiv != NULL)
     275            Xa->suiv->prec = Xa;
     276        }
     277        else {
     278          Xa->prec = X->prec;
     279          if (Xa->prec != NULL)
     280            Xa->prec->suiv = Xa;
     281        }
     282        Xa->inpivot = X->inpivot;
     283        if (X->inpivot != -1)        /* ecrase X dans pivot */
     284          pivot[X->inpivot] = Xa;
     285        Xa->inmodule = X->inmodule;
     286        if (X->inmodule != -1)        /* ecrase X dans pivot */
     287          module[X->inmodule] = Xa;
     288
     289        Xa->whereXa = 0;
     290        continue;
    291291      }
    292292
    293293      /* Maintenant on fait addpivot(X,Xa)
    294         noter que X et Xa sont non vides */
     294        noter que X et Xa sont non vides */
    295295
    296296      if (X->inpivot == -1) {
    297         if ((X->inmodule != -1)
    298             && (X->fin - X->debut < Xa->fin - Xa->debut)) {
    299           /* remplace X par Xa dans module */
    300           module[X->inmodule] = Xa;
    301           Xa->inmodule = X->inmodule;
    302           X->inmodule = -1;
    303           if (DEBUG)
    304             printf("Dans module %i-%i ecrase %i-%i\n",
    305                    1 + S[Xa->debut]->nom, 1 + S[Xa->fin]->nom,
    306                    1 + S[X->debut]->nom, 1 + S[X->fin]->nom);
    307         }
    308         else {
    309           if (X->inmodule == -1) {
    310             if (X->fin - X->debut < Xa->fin - Xa->debut)
    311               Z = Xa;
    312             else
    313               Z = X;
    314             /* ajoute Z (=max(X,Xa)) a module */
    315             module[(*imodule)] = Z;
    316             Z->inmodule = (*imodule);
    317             (*imodule)++;
    318             if (DEBUG)
    319               printf("module empile:%i-%i\n",
    320                      1 + S[Z->debut]->nom, 1 + S[Z->fin]->nom);
    321           }
    322         }
     297        if ((X->inmodule != -1)
     298            && (X->fin - X->debut < Xa->fin - Xa->debut)) {
     299          /* remplace X par Xa dans module */
     300          module[X->inmodule] = Xa;
     301          Xa->inmodule = X->inmodule;
     302          X->inmodule = -1;
     303          if (DEBUG)
     304            printf("Dans module %i-%i ecrase %i-%i\n",
     305                   1 + S[Xa->debut]->nom, 1 + S[Xa->fin]->nom,
     306                   1 + S[X->debut]->nom, 1 + S[X->fin]->nom);
     307        }
     308        else {
     309          if (X->inmodule == -1) {
     310            if (X->fin - X->debut < Xa->fin - Xa->debut)
     311              Z = Xa;
     312            else
     313              Z = X;
     314            /* ajoute Z (=max(X,Xa)) a module */
     315            module[(*imodule)] = Z;
     316            Z->inmodule = (*imodule);
     317            (*imodule)++;
     318            if (DEBUG)
     319              printf("module empile:%i-%i\n",
     320                     1 + S[Z->debut]->nom, 1 + S[Z->fin]->nom);
     321          }
     322        }
    323323      }
    324324
    325325      if (X->inpivot != -1)
    326         Z = Xa;
     326        Z = Xa;
    327327      else if (X->fin - X->debut < Xa->fin - Xa->debut)
    328         Z = X;
     328        Z = X;
    329329      else
    330         Z = Xa;
     330        Z = Xa;
    331331      /* on empile Z dans pivot */
    332332      pivot[(*ipivot)] = Z;
    333333      Z->inpivot = (*ipivot);
    334334      (*ipivot)++;
    335335      if (DEBUG)
    336         printf("pivot empile: %i-%i\n", 1 + S[Z->debut]->nom,
    337                1 + S[Z->fin]->nom);
     336        printf("pivot empile: %i-%i\n", 1 + S[Z->debut]->nom,
     337               1 + S[Z->fin]->nom);
    338338      X->whereXa = 0;
    339339      Xa->whereXa = 0;
    340340    }
     
    345345}
    346346
    347347sommet **algo1(graphe G)
    348      /* Entree: un graphe G 
    349         Sortie: une permutation factorisante de G,
    350         donnee sous la forme d'un tableau de structures Sommet ordonnees selon sigma.
    351         d'apres le travail de Habib/Paul/Viennot */
     348     /* Entree: un graphe G
     349        Sortie: une permutation factorisante de G,
     350        donnee sous la forme d'un tableau de structures Sommet ordonnees selon sigma.
     351        d'apres le travail de Habib/Paul/Viennot */
    352352{
    353353  int n; // nombre de sommets de G
    354354
    355   sclasse **pivot;              /*pile des pivots */
    356   int ipivot = 0;               /*indice sur la precedante */
     355  sclasse **pivot;                /*pile des pivots */
     356  int ipivot = 0;                /*indice sur la precedante */
    357357
    358     sclasse **module;           /*idem, modules */
     358    sclasse **module;                /*idem, modules */
    359359    int imodule = 0;
    360360
    361361    sclasse *singclasse;
    362362    /*invariant: toute classe avant singclasse dans la chaine */
    363363    /*a un seul element */
    364     int numclasse;              /* quand vaut n, on a fini!! */
     364    int numclasse;                /* quand vaut n, on a fini!! */
    365365
    366     sclasse *C1;                /*premiere classe, tete de la chaine */
    367     sclasse *Y;                 /*classe qui raffine */
    368     sclasse *X;                 /*classe raffinee */
    369     sclasse *Xa, *Xc;           /* morceaux de X */
    370     sommet *x;                  /* x in X */
    371     sommet *y;                  /* y in Y */
    372     sommet *centre;             /* le centre du raffinage actuel */
     366    sclasse *C1;                /*premiere classe, tete de la chaine */
     367    sclasse *Y;                        /*classe qui raffine */
     368    sclasse *X;                        /*classe raffinee */
     369    sclasse *Xa, *Xc;                /* morceaux de X */
     370    sommet *x;                        /* x in X */
     371    sommet *y;                        /* y in Y */
     372    sommet *centre;                /* le centre du raffinage actuel */
    373373
    374     sommet **S;                 /*la permutation factorisante ! */
     374    sommet **S;                        /*la permutation factorisante ! */
    375375
    376     int i, j;                   /*divers indices */
    377     sommet *scourant;           /* pour l'init */
    378     sadj *nextadj;              /*sommet adjacent suivant */
     376    int i, j;                        /*divers indices */
     377    sommet *scourant;                /* pour l'init */
     378    sadj *nextadj;                /*sommet adjacent suivant */
    379379    adj *nextadj2;              /* idem mais de type adj */
    380     info Inf;                   /* diverses info a passer a raffiner */
     380    info Inf;                        /* diverses info a passer a raffiner */
    381381
    382382    /* debut des initialisations */
    383383    n=G.n;
     
    393393    C1->debut = 0;
    394394    C1->fin = n - 1;
    395395    for (i = 0; i < n; i++) {
    396         /* initialisation des sommets */
    397         /* notre bebe est le sommet i dans M */
    398         scourant = (sommet *) fabmalloc(sizeof(struct Sommet));
    399         scourant->nom = i;
    400         scourant->place = i;    /* a ce point S=identite */
    401         scourant->adj = NULL; /* pas encore d'adjacence */
    402         scourant->classe = C1;
    403         S[i] = scourant;
     396        /* initialisation des sommets */
     397        /* notre bebe est le sommet i dans M */
     398        scourant = (sommet *) fabmalloc(sizeof(struct Sommet));
     399        scourant->nom = i;
     400        scourant->place = i;        /* a ce point S=identite */
     401        scourant->adj = NULL; /* pas encore d'adjacence */
     402        scourant->classe = C1;
     403        S[i] = scourant;
    404404    }
    405405    for (i = 0; i < n; i++)
    406406      {
    407         nextadj2 = G.G[i];
    408         while(nextadj2 != NULL)
    409           {
    410             j=nextadj2->s; //numero du sommet pointe
    411             if((j<0)||(j>=n))
    412               {
    413                 perror("Graphe invalide (numero de sommet erronne)!\n");
    414                 exit(1);
    415               }
    416             nextadj = (sadj *) fabmalloc(sizeof(struct Sadj));
    417             //un nouveau sadj
    418             nextadj->pointe = S[j];
    419             nextadj->suiv = S[i]->adj; //tete de liste
    420             if(nextadj->suiv!=NULL)
    421               nextadj->suiv->prec=nextadj;
    422             nextadj->prec=NULL;
    423             S[i]->adj = nextadj;        /*et le tour est joue */
    424             nextadj2=nextadj2->suiv;
    425           }
    426         }
     407        nextadj2 = G.G[i];
     408        while(nextadj2 != NULL)
     409          {
     410            j=nextadj2->s; //numero du sommet pointe
     411            if((j<0)||(j>=n))
     412              {
     413                perror("Graphe invalide (numero de sommet erronne)!\n");
     414                exit(1);
     415              }
     416            nextadj = (sadj *) fabmalloc(sizeof(struct Sadj));
     417            //un nouveau sadj
     418            nextadj->pointe = S[j];
     419            nextadj->suiv = S[i]->adj; //tete de liste
     420            if(nextadj->suiv!=NULL)
     421              nextadj->suiv->prec=nextadj;
     422            nextadj->prec=NULL;
     423            S[i]->adj = nextadj;        /*et le tour est joue */
     424            nextadj2=nextadj2->suiv;
     425          }
     426        }
    427427    /* NB: module et pivot sont vides */
    428428    Inf.pivot = pivot;
    429429    Inf.ipivot = &ipivot;
     
    434434    /* init terminnee */
    435435
    436436    while (1) {
    437         while (ipivot > 0 || imodule > 0) {
    438             while (ipivot > 0) {
    439                 /*cette boucle raffine selon tous les sommets
    440                    de la premiere classe dans pivot */
     437        while (ipivot > 0 || imodule > 0) {
     438            while (ipivot > 0) {
     439                /*cette boucle raffine selon tous les sommets
     440                   de la premiere classe dans pivot */
    441441
    442                 Y = pivot[ipivot - 1];
    443                 ipivot--;
    444                 Y->inpivot = -1;
     442                Y = pivot[ipivot - 1];
     443                ipivot--;
     444                Y->inpivot = -1;
    445445
    446                 for (i = Y->debut; i <= Y->fin; i++)
    447                     Raffiner(S, S[i], centre, &Inf);
     446                for (i = Y->debut; i <= Y->fin; i++)
     447                    Raffiner(S, S[i], centre, &Inf);
    448448
    449                 /* une optimisation de la fin de l'algo */
    450                 if (numclasse == n)
    451                     return (S);
    452             }
    453             /*maintenant pivot est vide, mais peut-etre pas module */
    454             if (imodule > 0) {
    455                 /* relance par un sommet (pas au pif...) */
    456                 /* de chaque module qui le represente */
    457                 Y = module[imodule - 1];
    458                 imodule--;
    459                 Y->inmodule = -1;
    460                 y = S[Y->debut];        /* le firstpivot sera toujours... */
    461                 Y->firstpivot = y;      /* le premier!! */
    462                 if (DEBUG)
    463                     printf("module-pivot %i-%i: sommet %i\n",
    464                            1 + S[Y->debut]->nom, 1 + S[Y->fin]->nom,
    465                            1 + y->nom);
    466                 Raffiner(S, y, centre, &Inf);
    467             }
    468         }
    469         /* a ce point, pivot et module sont vides...
    470            pas de pb! On va faire initpartition HERE */
    471         if (DEBUG)
    472             printf("\nInit Partition\n");
    473         /**** ajoute ici pour debbugger, mais moche!! */
    474         singclasse = S[0]->classe;
    475         while ((singclasse != NULL) &&
    476                (singclasse->debut == singclasse->fin))
    477           {
    478             singclasse = singclasse->suiv;
    479           }
    480         /* singclasse est la premiere classe
    481            non singlette, sauf si: */
    482         if (singclasse == NULL)
    483             /* on a n classes singlettes? ben c'est gagne! */
    484           {
    485             return (S);
    486           }
    487         if (singclasse == NULL && numclasse < n) {
    488             perror("c'est pas normal! Ca termine trop vite!\n");
    489             exit(1);
    490         }
     449                /* une optimisation de la fin de l'algo */
     450                if (numclasse == n)
     451                    return (S);
     452            }
     453            /*maintenant pivot est vide, mais peut-etre pas module */
     454            if (imodule > 0) {
     455                /* relance par un sommet (pas au pif...) */
     456                /* de chaque module qui le represente */
     457                Y = module[imodule - 1];
     458                imodule--;
     459                Y->inmodule = -1;
     460                y = S[Y->debut];        /* le firstpivot sera toujours... */
     461                Y->firstpivot = y;        /* le premier!! */
     462                if (DEBUG)
     463                    printf("module-pivot %i-%i: sommet %i\n",
     464                           1 + S[Y->debut]->nom, 1 + S[Y->fin]->nom,
     465                           1 + y->nom);
     466                Raffiner(S, y, centre, &Inf);
     467            }
     468        }
     469        /* a ce point, pivot et module sont vides...
     470           pas de pb! On va faire initpartition HERE */
     471        if (DEBUG)
     472            printf("\nInit Partition\n");
     473        /**** ajoute ici pour debbugger, mais moche!! */
     474        singclasse = S[0]->classe;
     475        while ((singclasse != NULL) &&
     476               (singclasse->debut == singclasse->fin))
     477          {
     478            singclasse = singclasse->suiv;
     479          }
     480        /* singclasse est la premiere classe
     481           non singlette, sauf si: */
     482        if (singclasse == NULL)
     483            /* on a n classes singlettes? ben c'est gagne! */
     484          {
     485            return (S);
     486          }
     487        if (singclasse == NULL && numclasse < n) {
     488            perror("c'est pas normal! Ca termine trop vite!\n");
     489            exit(1);
     490        }
    491491
    492         X = singclasse;
    493         x = X->firstpivot;
    494         if (x == NULL)
    495             x = S[X->debut];
    496         else                    /* remet firstpivot a NULL!! */
    497             X->firstpivot = NULL;
     492        X = singclasse;
     493        x = X->firstpivot;
     494        if (x == NULL)
     495            x = S[X->debut];
     496        else                        /* remet firstpivot a NULL!! */
     497            X->firstpivot = NULL;
    498498
    499         if (DEBUG)
    500             printf("Relance dans le module %i-%i avec le sommet %i\n",
    501                    1 + S[X->debut]->nom, 1 + S[X->fin]->nom, 1 + x->nom);
     499        if (DEBUG)
     500            printf("Relance dans le module %i-%i avec le sommet %i\n",
     501                   1 + S[X->debut]->nom, 1 + S[X->fin]->nom, 1 + x->nom);
    502502
    503         centre = x;             /*important! */
    504         /* astuce: on place {x} en tete de X
    505            ensuite, on raffine S selon x -> seule X est coupee
    506            il y a alors {x} X Xa
    507            -> on met {x} en queue de X et c'est bon!
    508            ainsi on a bien nonvoisins-x-voisons */
    509         Xc = nouvclasse(X->prec, X);
    510         numclasse++;
    511         x->classe = Xc;
    512         permute(S, x->place, X->debut);
    513         X->debut++;
    514         Xc->debut = x->place;
    515         Xc->fin = x->place;
    516         Raffiner(S, x, x, &Inf);
    517         /* X existe-il encore? */
    518         if (X->debut > X->fin)
    519             continue;
    520         /* echange de x et {x}. Init: -{x}-X- */
    521         Xc->suiv = X->suiv;
    522         if (X->suiv != NULL)
    523             X->suiv->prec = Xc;
    524         X->prec = Xc->prec;
    525         if (Xc->prec != NULL)
    526             Xc->prec->suiv = X;
    527         X->suiv = Xc;
    528         Xc->prec = X;
    529         permute(S, x->place, X->fin);
    530         Xc->debut = x->place;
    531         Xc->fin = x->place;
    532         X->debut--;
    533         X->fin--;
    534         //antibug?
    535         singclasse=X;
    536         /* now -X-{x}- */
    537         if (DEBUG)
    538             printS(S, n);
     503        centre = x;                /*important! */
     504        /* astuce: on place {x} en tete de X
     505           ensuite, on raffine S selon x -> seule X est coupee
     506           il y a alors {x} X Xa
     507           -> on met {x} en queue de X et c'est bon!
     508           ainsi on a bien nonvoisins-x-voisons */
     509        Xc = nouvclasse(X->prec, X);
     510        numclasse++;
     511        x->classe = Xc;
     512        permute(S, x->place, X->debut);
     513        X->debut++;
     514        Xc->debut = x->place;
     515        Xc->fin = x->place;
     516        Raffiner(S, x, x, &Inf);
     517        /* X existe-il encore? */
     518        if (X->debut > X->fin)
     519            continue;
     520        /* echange de x et {x}. Init: -{x}-X- */
     521        Xc->suiv = X->suiv;
     522        if (X->suiv != NULL)
     523            X->suiv->prec = Xc;
     524        X->prec = Xc->prec;
     525        if (Xc->prec != NULL)
     526            Xc->prec->suiv = X;
     527        X->suiv = Xc;
     528        Xc->prec = X;
     529        permute(S, x->place, X->fin);
     530        Xc->debut = x->place;
     531        Xc->fin = x->place;
     532        X->debut--;
     533        X->fin--;
     534        //antibug?
     535        singclasse=X;
     536        /* now -X-{x}- */
     537        if (DEBUG)
     538            printS(S, n);
    539539    }
    540540}
    541541
    542542/***************************************************************
    543 Etape intermediaire: trier toutes les listes d'adjacence 
     543Etape intermediaire: trier toutes les listes d'adjacence
    544544selon S. ce sont les listes de type sadj qui sont concernees
    545545***************************************************************/
    546546int Calculm(graphe G)
     
    552552    {
    553553      a=G.G[i];
    554554      while(a!=NULL)
    555         {
    556           a=a->suiv;
    557           r++;
    558         }
     555        {
     556          a=a->suiv;
     557          r++;
     558        }
    559559    }
    560   if(r%2!=0) 
    561     { 
     560  if(r%2!=0)
     561    {
    562562      perror("Erreur: nombre impaire d'arrete, graphe non-oriente??\n");
    563563      exit(1);
    564564    }
     
    584584    {
    585585      a=S[i]->adj;
    586586      while(a!=NULL)
    587         {
    588           tab1[i]++;
    589           a=a->suiv;
    590         }
     587        {
     588          tab1[i]++;
     589          a=a->suiv;
     590        }
    591591    }
    592592  //deuxieme passe: frequences cumulees a rebours
    593593  // (car les listes d'adjacences se construisent a l'envers
     
    595595  //for(i=n-1;i>0;i--)
    596596  //  tab1[i-1]+=tab1[i];
    597597
    598   //deuxieme passe: frequences cumulees 
     598  //deuxieme passe: frequences cumulees
    599599  for(i=1;i<n;i++)
    600600    tab1[i]+=tab1[i-1];
    601  
     601
    602602  //troisieme passe: liste double
    603603  for(i=0; i<n; i++)
    604604    {
    605605      a=S[i]->adj;
    606606      while(a!=NULL)
    607         {
    608           /* cree un nouveau record */
    609           c=(clef2 *)fabmalloc(sizeof(struct Clef2));
    610           c->i=i;
    611           c->nom=a->pointe->nom;
    612           c->place=a->pointe->place;
    613           /* le place bien dans tab2 */
    614           tab1[c->place]--;
    615           tab2[tab1[c->place]]=c;
    616           /*et on continue */
    617           a=a->suiv; 
    618         }
     607        {
     608          /* cree un nouveau record */
     609          c=(clef2 *)fabmalloc(sizeof(struct Clef2));
     610          c->i=i;
     611          c->nom=a->pointe->nom;
     612          c->place=a->pointe->place;
     613          /* le place bien dans tab2 */
     614          tab1[c->place]--;
     615          tab2[tab1[c->place]]=c;
     616          /*et on continue */
     617          a=a->suiv;
     618        }
    619619    }
    620620
    621621  //quatrieme passe: detruit les vielles listes d'adjacence
     
    623623    {
    624624      a=S[i]->adj;
    625625      while(a!=NULL)
    626         {
    627           atmp=a->suiv;
    628           free(a);
    629           a=atmp;
    630         }
     626        {
     627          atmp=a->suiv;
     628          free(a);
     629          a=atmp;
     630        }
    631631      S[i]->adj=NULL;
    632632    }
    633633
     
    639639      a->pointe=S[c->i];
    640640      a->suiv=S[c->place]->adj; //insere en tete
    641641      if(a->suiv!=NULL)
    642         a->suiv->prec=a;
     642        a->suiv->prec=a;
    643643      a->prec=NULL;
    644644      S[c->place]->adj=a;
    645645      //nettoie
    646       free(c); 
     646      free(c);
    647647   }
    648648  free(tab1);
    649649  free(tab2);
     
    676676    /*   de "abscence de separateur" (-1). De plus, on fera des min et des */
    677677    /*   max sur les bords */
    678678    nn->bg = n + 2;
    679     nn->bd = -2;                /* idem */
     679    nn->bd = -2;                /* idem */
    680680
    681681    nn->fils = NULL;
    682682    nn->lastfils = NULL;
     
    693693    nf->pointe = nfils;
    694694    nf->suiv = NULL;
    695695    if (pere->fils == NULL)
    696         pere->fils = nf;        /* on cree le premier fils */
     696        pere->fils = nf;        /* on cree le premier fils */
    697697    else
    698         pere->lastfils->suiv = nf;      /* on ajoute nf a la chaine */
     698        pere->lastfils->suiv = nf;        /* on ajoute nf a la chaine */
    699699    pere->lastfils = nf;
    700     nfils->pere = pere;         /* normalement: redondant,mais... */
     700    nfils->pere = pere;                /* normalement: redondant,mais... */
    701701    nfils->fpere = nf;
    702702}
    703703
     
    711711    /* met a jour la liste des peres */
    712712    f = artefact->fils;
    713713    while (f != NULL) {
    714         f->pointe->pere = pere; /*avant c'etait ancien... */
    715         /* f->pointe->fpere est inchange */
    716         f = f->suiv;
     714        f->pointe->pere = pere;        /*avant c'etait ancien... */
     715        /* f->pointe->fpere est inchange */
     716        f = f->suiv;
    717717    }
    718718    /* greffe la liste */
    719719    greffe = artefact->fpere;
    720720    artefact->lastfils->suiv = greffe->suiv;
    721721    greffe->pointe = artefact->fils->pointe;
    722722    greffe->suiv = artefact->fils->suiv;
    723     artefact->fils->pointe->fpere = greffe;     /*artefact->fils a disparu */
     723    artefact->fils->pointe->fpere = greffe;        /*artefact->fils a disparu */
    724724    if (pere->lastfils == greffe)
    725         pere->lastfils = artefact->lastfils;
     725        pere->lastfils = artefact->lastfils;
    726726}
    727727
    728728void
     
    730730{
    731731    /* extrait la liste [premier...dernier] des fils de l'ancien noeud,
    732732       et en fait la liste des fils du nouveau noeud */
    733     fils *nf;                   /* il faut une structure fils de plus */
    734     fils *f;                    /*indice de mise a jour */
     733    fils *nf;                        /* il faut une structure fils de plus */
     734    fils *f;                        /*indice de mise a jour */
    735735    nf = (fils *) fabmalloc(sizeof(fils));
    736736    nf->pointe = premier->pointe;
    737737    nf->suiv = premier->suiv;
     
    741741    nouveau->lastfils = dernier;
    742742    nouveau->bg = premier->pointe->bg;
    743743    nouveau->bd = dernier->pointe->bd;
    744     nouveau->ps = premier->pointe->bg;  /* nouveau est suppose etre un */
    745     nouveau->ds = dernier->pointe->bd;  /* module, donc bords=separateurs! */
     744    nouveau->ps = premier->pointe->bg;        /* nouveau est suppose etre un */
     745    nouveau->ds = dernier->pointe->bd;        /* module, donc bords=separateurs! */
    746746    if (ancien->lastfils == dernier)
    747         ancien->lastfils = premier;
     747        ancien->lastfils = premier;
    748748    /* ecrase l'ancier premier */
    749749    nouveau->fpere = premier;
    750750    premier->pointe = nouveau;
     
    754754    /* met a jour la liste des peres */
    755755    f = nf;
    756756    while (f != dernier->suiv) {
    757         f->pointe->pere = nouveau;      /*avant c'etait ancien... */
    758         f->pointe->fpere = premier;
    759         f = f->suiv;
     757        f->pointe->pere = nouveau;        /*avant c'etait ancien... */
     758        f->pointe->fpere = premier;
     759        f = f->suiv;
    760760    }
    761761}
    762762
     
    769769    ffils = N->fils;
    770770
    771771    for (i = 0; i < level - 1; i++)
    772         printf("  |");
     772        printf("  |");
    773773    if (N->pere == NULL)
    774         printf(" ");
     774        printf(" ");
    775775    else
    776         printf("  +-");
     776        printf("  +-");
    777777    switch (N->type) {
    778778    case UNKN:
    779         printf("Noeud\n");
    780         break;
     779        printf("Noeud\n");
     780        break;
    781781    case MODULE:
    782         printf("Module\n");
    783         break;
     782        printf("Module\n");
     783        break;
    784784    case ARTEFACT:
    785         printf("Artefact\n");
    786         break;
     785        printf("Artefact\n");
     786        break;
    787787    case SERIE:
    788         printf("Serie \n");
    789         break;
     788        printf("Serie \n");
     789        break;
    790790    case PARALLELE:
    791         printf("Parallele \n");
    792         break;
     791        printf("Parallele \n");
     792        break;
    793793    case PREMIER:
    794         printf("Premier \n");
    795         break;
     794        printf("Premier \n");
     795        break;
    796796    }
    797797
    798798    do {
    799         nfils = ffils->pointe;
    800         if (nfils->type == FEUILLE) {
    801             for (i = 0; i < level; i++)
    802                 printf("  |");
    803             printf("  +--");
    804             printf("%i\n", 1 + nfils->nom);
    805         }
    806         else {
    807             printnoeud(nfils, level + 1);
    808         }
    809         ffils = ffils->suiv;
     799        nfils = ffils->pointe;
     800        if (nfils->type == FEUILLE) {
     801            for (i = 0; i < level; i++)
     802                printf("  |");
     803            printf("  +--");
     804            printf("%i\n", 1 + nfils->nom);
     805        }
     806        else {
     807            printnoeud(nfils, level + 1);
     808        }
     809        ffils = ffils->suiv;
    810810    }
    811811    while (ffils != NULL);
    812812}
     
    824824*/
    825825    /* debug: S n'est utilise que pour mettre vrainom a jour */
    826826    int n; //nombre de sommets du graphe
    827     int *ouvrantes;             /* tableau du nombre de parentheses ouvrantes */
     827    int *ouvrantes;                /* tableau du nombre de parentheses ouvrantes */
    828828    /* ouvrante[i]=3 ssi i-1(((i  */
    829829    /* ouvrante [0]=3: (((0 */
    830830
    831     int *fermantes;             /* idem fermantes[i]=2 ssi i)))i+1
    832                                    fermante [n-1]=2 ssi n))) */
    833     int *ps;                    /* ps[i]=premier separateur de (i,i+1) */
     831    int *fermantes;                /* idem fermantes[i]=2 ssi i)))i+1
     832                                   fermante [n-1]=2 ssi n))) */
     833    int *ps;                        /* ps[i]=premier separateur de (i,i+1) */
    834834    int *ds;
    835835
    836     int i, j;                   /*indices de paires ou de sommets */
     836    int i, j;                        /*indices de paires ou de sommets */
    837837
    838838    sadj *a1, *a2;                /* parcours de liste d'adjacence */
    839839
    840     noeud *racine;              /*racine du pseudocoardre */
    841     noeud *courant, *nouveau;   /* noeud courant du pseudocoarbre */
    842     noeud **pileinterne;        /* pile des modules pour les passes 3,5,5 */
    843     int indicepileinterne = 0;  /*pointeur dans cette pile */
    844     int taillepileinterne;      /* taille de la pile apres la 2eme passe */
     840    noeud *racine;                /*racine du pseudocoardre */
     841    noeud *courant, *nouveau;        /* noeud courant du pseudocoarbre */
     842    noeud **pileinterne;        /* pile des modules pour les passes 3,5,5 */
     843    int indicepileinterne = 0;        /*pointeur dans cette pile */
     844    int taillepileinterne;        /* taille de la pile apres la 2eme passe */
    845845
    846846    int *adjii;                 /* adjii[i]=1 ssi S[i] et S[i+1] sont */
    847                                    /*                      adjacents */
     847                                   /*                             adjacents */
    848848    /*PROPHASE : initialisations */
    849849    n=G.n;
    850850    ouvrantes = (int *) fabmalloc(n * sizeof(int));
     
    863863    /* remplit adjii qui dit si S[i] adjacent a S[i+1] */
    864864    for(i=0; i<n-1; i++)
    865865      {
    866         a1=S[i]->adj;
    867         while((a1!=NULL)&&(a1->pointe->place != i+1))
    868            a1=a1->suiv;
    869         if( a1 == NULL)
    870            adjii[i]=0;
    871         else // a1->pointe->place==i+1, donc i adj i+1
    872            adjii[i]=1;
     866        a1=S[i]->adj;
     867        while((a1!=NULL)&&(a1->pointe->place != i+1))
     868           a1=a1->suiv;
     869        if( a1 == NULL)
     870           adjii[i]=0;
     871        else // a1->pointe->place==i+1, donc i adj i+1
     872           adjii[i]=1;
    873873      }
    874874    adjii[n-1]=0; //perfectionnisme
    875875
     
    880880       complexite: O(n^2) */
    881881
    882882    ouvrantes[0] = 1;
    883     fermantes[n - 1] = 1;       /* parentheses des bords */
     883    fermantes[n - 1] = 1;        /* parentheses des bords */
    884884
    885885    for (i = 0; i < n - 1; i++) {
    886886      /*recherche de ps(i,i+1) */
    887887      a1=S[i]->adj;
    888888      a2=S[i+1]->adj;
    889889      while((a1!=NULL) && (a2!=NULL) && (a1->pointe->place<i) &&
    890             (a2->pointe->place<i) && (a1->pointe->place == a2->pointe->place))
    891         {
    892           a1=a1->suiv;
    893           a2=a2->suiv;
    894         }
     890            (a2->pointe->place<i) && (a1->pointe->place == a2->pointe->place))
     891        {
     892          a1=a1->suiv;
     893          a2=a2->suiv;
     894        }
    895895
    896896      //arbre de decision complique pour trouver le premier separateur!
    897       if( ((a1==NULL) && (a2==NULL)) 
    898           ||((a1==NULL) &&(a2->pointe->place >= i))
    899           ||((a2==NULL) && (a1->pointe->place >= i))
    900           ||((a1!=NULL) && (a2!=NULL) && (a1->pointe->place >= i) && (a2->pointe->place >= i)))
    901         //pas de separateur
    902         ps[i]=i+1;
     897      if( ((a1==NULL) && (a2==NULL))
     898          ||((a1==NULL) &&(a2->pointe->place >= i))
     899          ||((a2==NULL) && (a1->pointe->place >= i))
     900          ||((a1!=NULL) && (a2!=NULL) && (a1->pointe->place >= i) && (a2->pointe->place >= i)))
     901        //pas de separateur
     902        ps[i]=i+1;
    903903      else
    904         {
    905           if((a1==NULL) || (a1->pointe->place >= i))
    906             ps[i]=a2->pointe->place;
    907           else if((a2==NULL) || (a2->pointe->place >= i))
    908             ps[i]=a1->pointe->place;
    909           else
    910             {
    911               if((a1->suiv!=NULL)&&(a1->suiv->pointe->place == a2->pointe->place))
    912                 ps[i]=a1->pointe->place;
    913               else if ((a2->suiv!=NULL)&&(a2->suiv->pointe->place == a1->pointe->place))
    914                 ps[i]=a2->pointe->place;
    915               else
    916                 ps[i]=min(a1->pointe->place , a2->pointe->place);
    917             }
    918           ouvrantes[ps[i]]++;   /* marque la fracture gauche, if any */
    919           fermantes[i]++;
    920         }
     904        {
     905          if((a1==NULL) || (a1->pointe->place >= i))
     906            ps[i]=a2->pointe->place;
     907          else if((a2==NULL) || (a2->pointe->place >= i))
     908            ps[i]=a1->pointe->place;
     909          else
     910            {
     911              if((a1->suiv!=NULL)&&(a1->suiv->pointe->place == a2->pointe->place))
     912                ps[i]=a1->pointe->place;
     913              else if ((a2->suiv!=NULL)&&(a2->suiv->pointe->place == a1->pointe->place))
     914                ps[i]=a2->pointe->place;
     915              else
     916                ps[i]=min(a1->pointe->place , a2->pointe->place);
     917            }
     918          ouvrantes[ps[i]]++;        /* marque la fracture gauche, if any */
     919          fermantes[i]++;
     920        }
    921921      if (DEBUG)
    922         printf("ps(%i,%i)=%i\n", i , i+1, ps[i]);
     922        printf("ps(%i,%i)=%i\n", i , i+1, ps[i]);
    923923
    924         /*recherche de ds(i,i+1)
    925           plus penible encore!*/
     924        /*recherche de ds(i,i+1)
     925          plus penible encore!*/
    926926      a1=S[i]->adj;
    927927      if(a1!=NULL) // se place en queue de liste.
    928         while(a1->suiv!=NULL)
    929           a1=a1->suiv;
     928        while(a1->suiv!=NULL)
     929          a1=a1->suiv;
    930930      a2=S[i+1]->adj;
    931931      if(a2!=NULL)
    932         while(a2->suiv!=NULL)
    933           a2=a2->suiv;
     932        while(a2->suiv!=NULL)
     933          a2=a2->suiv;
    934934      while((a1!=NULL) && (a2!=NULL) && (a1->pointe->place > i+1) &&
    935             (a2->pointe->place > i+1) && (a1->pointe->place == a2->pointe->place))
    936         {
    937           a1=a1->prec;
    938           a2=a2->prec;
    939         }
    940       if( ((a1==NULL) && (a2==NULL)) 
    941           ||((a1==NULL) && (a2->pointe->place <= i+1))
    942           ||((a2==NULL) && (a1->pointe->place <= i+1))
    943           ||((a1!=NULL) && (a2!=NULL) && (a1->pointe->place <= i+1) && (a2->pointe->place <= i+1)))
    944         //pas de separateur
    945         ds[i]=i+1;
     935            (a2->pointe->place > i+1) && (a1->pointe->place == a2->pointe->place))
     936        {
     937          a1=a1->prec;
     938          a2=a2->prec;
     939        }
     940      if( ((a1==NULL) && (a2==NULL))
     941          ||((a1==NULL) && (a2->pointe->place <= i+1))
     942          ||((a2==NULL) && (a1->pointe->place <= i+1))
     943          ||((a1!=NULL) && (a2!=NULL) && (a1->pointe->place <= i+1) && (a2->pointe->place <= i+1)))
     944        //pas de separateur
     945        ds[i]=i+1;
    946946      else
    947         {
    948           if((a1==NULL) || (a1->pointe->place <= i+1))
    949             ds[i]=a2->pointe->place;
    950           else if((a2==NULL) || (a2->pointe->place <= i+1))
    951             ds[i]=a1->pointe->place;
    952           else
    953             {
    954               if((a1->prec!=NULL)&&(a1->prec->pointe->place == a2->pointe->place))
    955                 ds[i]=a1->pointe->place;
    956               else if((a2->prec!=NULL)&&(a2->prec->pointe->place == a1->pointe->place))
    957                 ds[i]=a2->pointe->place;
    958               else
    959                 ds[i]=max(a1->pointe->place , a2->pointe->place);
    960             }
    961        
    962        
    963           //ds[i] = j;
    964         ouvrantes[i + 1]++;     /* marque la fracture gauche, if any */
    965         fermantes[ds[i]]++;     /* attention aux decalages d'indices */
    966         }
     947        {
     948          if((a1==NULL) || (a1->pointe->place <= i+1))
     949            ds[i]=a2->pointe->place;
     950          else if((a2==NULL) || (a2->pointe->place <= i+1))
     951            ds[i]=a1->pointe->place;
     952          else
     953            {
     954              if((a1->prec!=NULL)&&(a1->prec->pointe->place == a2->pointe->place))
     955                ds[i]=a1->pointe->place;
     956              else if((a2->prec!=NULL)&&(a2->prec->pointe->place == a1->pointe->place))
     957                ds[i]=a2->pointe->place;
     958              else
     959                ds[i]=max(a1->pointe->place , a2->pointe->place);
     960            }
     961
     962
     963          //ds[i] = j;
     964        ouvrantes[i + 1]++;        /* marque la fracture gauche, if any */
     965        fermantes[ds[i]]++;        /* attention aux decalages d'indices */
     966        }
    967967      if (DEBUG)
    968         printf("ds(%i,%i)=%i\n", i,i+1,ds[i]);
    969       //S[i]->nom + 1,         S[i + 1]->nom + 1, S[ds[i]]->nom + 1);
     968        printf("ds(%i,%i)=%i\n", i,i+1,ds[i]);
     969      //S[i]->nom + 1,               S[i + 1]->nom + 1, S[ds[i]]->nom + 1);
    970970    }
    971  
     971
    972972    /*DEUXIEME PASSE: construction du pseudocoarbre */
    973973
    974974    racine = nouvnoeud(UNKN, NULL, -1, n);
    975975    courant = racine;
    976976    for (i = 0; i < n; i++) {
    977         /*1: on lit des parentheses ouvrantes: descentes */
    978         for (j = 0; j < ouvrantes[i]; j++) {
    979             /*Descente vers un nouveau noeud */
    980             nouveau = nouvnoeud(UNKN, courant, -1, n);
    981             ajoutfils(courant, nouveau);        /*on l'ajoute... */
    982             courant = nouveau;  /* et on descent */
    983             if (DEBUG)
    984                 printf("(");
    985         }
    986         /* 2: on lit le sommet: feuille */
    987         nouveau = nouvnoeud(FEUILLE, courant, i, n);
    988         ajoutfils(courant, nouveau);    /*on ajoute le bebe... */
    989         (*nouveau).ps = i;
    990         (*nouveau).ds = i;
    991         (*nouveau).bg = i;
    992         (*nouveau).bd = i;
    993         nouveau->nom = S[i]->nom;
    994         /* et pourquoi pas i ? Je m'embrouille... */
    995         if (DEBUG)
    996             printf(" %i ", S[i]->nom + 1);
     977        /*1: on lit des parentheses ouvrantes: descentes */
     978        for (j = 0; j < ouvrantes[i]; j++) {
     979            /*Descente vers un nouveau noeud */
     980            nouveau = nouvnoeud(UNKN, courant, -1, n);
     981            ajoutfils(courant, nouveau);        /*on l'ajoute... */
     982            courant = nouveau;        /* et on descent */
     983            if (DEBUG)
     984                printf("(");
     985        }
     986        /* 2: on lit le sommet: feuille */
     987        nouveau = nouvnoeud(FEUILLE, courant, i, n);
     988        ajoutfils(courant, nouveau);        /*on ajoute le bebe... */
     989        (*nouveau).ps = i;
     990        (*nouveau).ds = i;
     991        (*nouveau).bg = i;
     992        (*nouveau).bd = i;
     993        nouveau->nom = S[i]->nom;
     994        /* et pourquoi pas i ? Je m'embrouille... */
     995        if (DEBUG)
     996            printf(" %i ", S[i]->nom + 1);
    997997
    998         /*3: on lit des parentheses fermantes: remontees */
    999         for (j = 0; j < fermantes[i]; j++) {
    1000             /*ASTUCE: ici on va en profiter pour supprimer les
    1001                noeuds a un fils, afin d'economiser une passe */
    1002             if (courant->fils == courant->lastfils) {   /*just one son */
    1003                 courant->pere->lastfils->pointe = courant->fils->pointe;
    1004                 courant->fils->pointe->pere = courant->pere;
    1005                 courant->fils->pointe->fpere = courant->pere->lastfils;
    1006                 /* explication: le dernier fils de courant.pere est
    1007                    actuellement courant himself. Il est donc pointe par
    1008                    courant.pere.lastfils.pointe. Il suffit de changer ce
    1009                    pointeur pour qu'il pointe maintenant non plus sur courant,
    1010                    mais sur l'unique fils de courant: courant.fils.pointe.
    1011                    Ouf! */
    1012                 /* NB: courant est maintenant deconnecte de l'arbre.
    1013                    on pourrait faire un free() mais bon... */
    1014             }
    1015             else {
    1016                 /*on empile ce noeud interne.
    1017                    L'ordre est celui de la postvisite d'un DFS */
    1018                 pileinterne[indicepileinterne] = courant;
    1019                 indicepileinterne++;
    1020             }
    1021             /* et dans tous les cas: on remonte! */
    1022             courant = courant->pere;
    1023             if (DEBUG)
    1024                 printf(")");
    1025         }
     998        /*3: on lit des parentheses fermantes: remontees */
     999        for (j = 0; j < fermantes[i]; j++) {
     1000            /*ASTUCE: ici on va en profiter pour supprimer les
     1001               noeuds a un fils, afin d'economiser une passe */
     1002            if (courant->fils == courant->lastfils) {        /*just one son */
     1003                courant->pere->lastfils->pointe = courant->fils->pointe;
     1004                courant->fils->pointe->pere = courant->pere;
     1005                courant->fils->pointe->fpere = courant->pere->lastfils;
     1006                /* explication: le dernier fils de courant.pere est
     1007                   actuellement courant himself. Il est donc pointe par
     1008                   courant.pere.lastfils.pointe. Il suffit de changer ce
     1009                   pointeur pour qu'il pointe maintenant non plus sur courant,
     1010                   mais sur l'unique fils de courant: courant.fils.pointe.
     1011                   Ouf! */
     1012                /* NB: courant est maintenant deconnecte de l'arbre.
     1013                   on pourrait faire un free() mais bon... */
     1014            }
     1015            else {
     1016                /*on empile ce noeud interne.
     1017                   L'ordre est celui de la postvisite d'un DFS */
     1018                pileinterne[indicepileinterne] = courant;
     1019                indicepileinterne++;
     1020            }
     1021            /* et dans tous les cas: on remonte! */
     1022            courant = courant->pere;
     1023            if (DEBUG)
     1024                printf(")");
     1025        }
    10261026    }
    10271027    if (DEBUG)
    1028         printf("\n");
     1028        printf("\n");
    10291029
    10301030    /* on enleve un ptit defaut */
    10311031    racine = racine->fils->pointe;
    10321032    racine->pere = NULL;
    10331033    racine->fpere = NULL;
    10341034    if (DEBUG) {
    1035         printf("Arbre apres la deuxieme passe:\n");
    1036         printnoeud(racine, 0);
     1035        printf("Arbre apres la deuxieme passe:\n");
     1036        printnoeud(racine, 0);
    10371037    }
    10381038
    10391039    /*TROISIEME PASSE */
     
    10461046
    10471047    taillepileinterne = indicepileinterne;
    10481048    for (indicepileinterne = 0; indicepileinterne < taillepileinterne;
    1049         indicepileinterne++) {
    1050         noeud *scanne;
    1051         fils *nextfils;
    1052         noeud *succourant;
    1053         /* scanne est le noeud (pere) que l'on regarde;
    1054            nextfils parcours sa liste de fils;
    1055            courant est le fils actuellement examine et
    1056            succourant=succ(courant) */
    1057         noeud *debutnoeud;
    1058         fils *debutfils;
    1059         /*deux variables utilise pour la recherche de jumeaux:
    1060            debut du bloc max */
     1049        indicepileinterne++) {
     1050        noeud *scanne;
     1051        fils *nextfils;
     1052        noeud *succourant;
     1053        /* scanne est le noeud (pere) que l'on regarde;
     1054           nextfils parcours sa liste de fils;
     1055           courant est le fils actuellement examine et
     1056           succourant=succ(courant) */
     1057        noeud *debutnoeud;
     1058        fils *debutfils;
     1059        /*deux variables utilise pour la recherche de jumeaux:
     1060           debut du bloc max */
    10611061
    1062         scanne = pileinterne[indicepileinterne];        /*he oui, la pile */
    1063         nextfils = scanne->fils;        /*on commence au premier fils */
    1064         do {
    1065             /*la boucle chiante... cf mon memoire de DEA */
    1066             courant = nextfils->pointe;
    1067             /* bords */
    1068             scanne->bg = min(scanne->bg, courant->bg);
    1069             scanne->bd = max(scanne->bd, courant->bd);
    1070             /*separateurs */
    1071             scanne->ps = min(scanne->ps, courant->ps);
    1072             if (scanne->fils->pointe != courant)
    1073                 /*ce n'est pas le premier fils */
    1074                 scanne->ps = min(scanne->ps, ps[(courant->bg) - 1]);
    1075             scanne->ds = max(scanne->ds, courant->ds);
    1076             if (scanne->lastfils->pointe != courant)
    1077                 /*ce n'est pas le dernier fils */
    1078                 scanne->ds = max(scanne->ds, ds[courant->bd]);
     1062        scanne = pileinterne[indicepileinterne];        /*he oui, la pile */
     1063        nextfils = scanne->fils;        /*on commence au premier fils */
     1064        do {
     1065            /*la boucle chiante... cf mon memoire de DEA */
     1066            courant = nextfils->pointe;
     1067            /* bords */
     1068            scanne->bg = min(scanne->bg, courant->bg);
     1069            scanne->bd = max(scanne->bd, courant->bd);
     1070            /*separateurs */
     1071            scanne->ps = min(scanne->ps, courant->ps);
     1072            if (scanne->fils->pointe != courant)
     1073                /*ce n'est pas le premier fils */
     1074                scanne->ps = min(scanne->ps, ps[(courant->bg) - 1]);
     1075            scanne->ds = max(scanne->ds, courant->ds);
     1076            if (scanne->lastfils->pointe != courant)
     1077                /*ce n'est pas le dernier fils */
     1078                scanne->ds = max(scanne->ds, ds[courant->bd]);
    10791079
    1080             nextfils = nextfils->suiv;
    1081         }
    1082         while (nextfils != NULL);
     1080            nextfils = nextfils->suiv;
     1081        }
     1082        while (nextfils != NULL);
    10831083
    10841084
    1085         if (DEBUG)
    1086             printf("noeud %i-%i: ps=%i ds=%i", 1 + scanne->bg,
    1087                    1 + scanne->bd, 1 + scanne->ps, 1 + scanne->ds);
     1085        if (DEBUG)
     1086            printf("noeud %i-%i: ps=%i ds=%i", 1 + scanne->bg,
     1087                   1 + scanne->bd, 1 + scanne->ps, 1 + scanne->ds);
    10881088
    1089         /* maintenant le test tout simple pour savoir si on a un module: */
    1090         if (((scanne->ps) == (scanne->bg))
    1091             && ((scanne->ds) == (scanne->bd))) {
    1092             /*on a un module */
    1093             scanne->type = MODULE;
    1094             if (DEBUG)
    1095                 printf(" Module.\n");
    1096         }
    1097         else {
    1098             scanne->type = ARTEFACT;
    1099             if (DEBUG)
    1100                 printf(" artefact.\n");
    1101         }
     1089        /* maintenant le test tout simple pour savoir si on a un module: */
     1090        if (((scanne->ps) == (scanne->bg))
     1091            && ((scanne->ds) == (scanne->bd))) {
     1092            /*on a un module */
     1093            scanne->type = MODULE;
     1094            if (DEBUG)
     1095                printf(" Module.\n");
     1096        }
     1097        else {
     1098            scanne->type = ARTEFACT;
     1099            if (DEBUG)
     1100                printf(" artefact.\n");
     1101        }
    11021102    }
    11031103
    11041104    if (DEBUG) {
    1105         printf("Arbre apres la troisieme passe:\n");
    1106         printnoeud(racine, 0);
     1105        printf("Arbre apres la troisieme passe:\n");
     1106        printnoeud(racine, 0);
    11071107    }
    11081108
    11091109    /* QUATRIEME PASSE */
     
    11111111       ca se fait de bas en haut grace a pileinterne (toujours elle!) */
    11121112
    11131113    for (indicepileinterne = 0; indicepileinterne < taillepileinterne;
    1114         indicepileinterne++) {
    1115         noeud *scanne;
    1116         scanne = pileinterne[indicepileinterne];        /*he oui, la pile */
    1117         if (scanne->type == ARTEFACT) {
    1118             /*attention! La pile peut contenir des noeuds deconnectes */
    1119             fusionne(scanne->pere, scanne);
    1120             if (DEBUG)
    1121                 printf("Artefact elimine: %i-%i\n", 1 + scanne->bg,
    1122                        1 + scanne->bd);
    1123         }
     1114        indicepileinterne++) {
     1115        noeud *scanne;
     1116        scanne = pileinterne[indicepileinterne];        /*he oui, la pile */
     1117        if (scanne->type == ARTEFACT) {
     1118            /*attention! La pile peut contenir des noeuds deconnectes */
     1119            fusionne(scanne->pere, scanne);
     1120            if (DEBUG)
     1121                printf("Artefact elimine: %i-%i\n", 1 + scanne->bg,
     1122                       1 + scanne->bd);
     1123        }
    11241124    }
    11251125    if (DEBUG) {
    1126         printf("Arbre apres la quatrieme passe:\n");
    1127         printnoeud(racine, 0);
     1126        printf("Arbre apres la quatrieme passe:\n");
     1127        printnoeud(racine, 0);
    11281128    }
    11291129
    11301130    /* CINQIEME ET DERNIERE PASSE */
    11311131    /* on va typer les noeuds et extraire les fusions.
    11321132       comment on fait? Ben, la pile.... */
    11331133    for (indicepileinterne = 0; indicepileinterne < taillepileinterne;
    1134         indicepileinterne++) {
    1135         noeud *scanne;
    1136         fils *nextfils;
    1137         noeud *succourant;
    1138         /* scanne est le noeud (pere) que l'on regarde;
    1139            nextfils parcours sa liste de fils;
    1140            courant est le fils actuellement examine et
    1141            succourant=succ(courant) */
    1142         noeud *debutnoeud;
    1143         fils *debutfils;
    1144         /*deux variables utilise pour la recherche de jumeaux:
    1145            debut du bloc max */
     1134        indicepileinterne++) {
     1135        noeud *scanne;
     1136        fils *nextfils;
     1137        noeud *succourant;
     1138        /* scanne est le noeud (pere) que l'on regarde;
     1139           nextfils parcours sa liste de fils;
     1140           courant est le fils actuellement examine et
     1141           succourant=succ(courant) */
     1142        noeud *debutnoeud;
     1143        fils *debutfils;
     1144        /*deux variables utilise pour la recherche de jumeaux:
     1145           debut du bloc max */
    11461146
    1147         scanne = pileinterne[indicepileinterne];
    1148         if (scanne->type != MODULE)
    1149             continue;           /* je traite que les modules */
     1147        scanne = pileinterne[indicepileinterne];
     1148        if (scanne->type != MODULE)
     1149            continue;                /* je traite que les modules */
    11501150
    1151         nextfils = scanne->fils;        /*on commence au premier fils */
    1152         while (1) {
    1153             courant = nextfils->pointe;
    1154             succourant = nextfils->suiv->pointe;
    1155             if (ps[courant->bd] >= courant->bg
    1156                 && ds[courant->bd] <= succourant->bd) {
    1157                 /*Des jumeaux!! ->serie ou parallele! */
    1158                 /* on va determiner le bloc max de jumeaux consecutifs */
    1159                 debutnoeud = courant;
    1160                 debutfils = nextfils;
    1161                 while (ps[courant->bd] >= courant->bg &&
    1162                        ds[courant->bd] <= succourant->bd &&
    1163                        nextfils->suiv != NULL) {
    1164                     nextfils = nextfils->suiv;
    1165                     courant = nextfils->pointe;
    1166                     if (nextfils->suiv == NULL)
    1167                         break;
    1168                     succourant = nextfils->suiv->pointe;
    1169                 }
    1170                 /*maintenant on connait la taille du bloc: il va de
    1171                    debutnoeud a courant inclus,
    1172                    et dans la liste des fils de scanne,
    1173                    il s'etend de debutfils a nextfils inclus.
    1174                    On extrait cette liste pour en faire les fils d'un
    1175                    nouveau noeud... sauf si pas la peine! */
    1176                 if (debutfils == scanne->fils
    1177                     && nextfils == scanne->lastfils) {
    1178                     /* le noeud scanne etait serie ou parallele */
    1179                   if ( adjii[debutnoeud->bd] !=0)
    1180                         scanne->type = SERIE;
    1181                     else
    1182                         scanne->type = PARALLELE;
    1183                 }
    1184                 else {
    1185                   if ( adjii[debutnoeud->bd]!=0)
    1186                     /*seule cette ligne distingue G de G' !! */
    1187                     {
    1188                         nouveau = nouvnoeud(SERIE, scanne, -1, n);
    1189                         if (DEBUG)
    1190                             printf("Module serie extrait: %i-%i\n",
    1191                                    1 + debutnoeud->bg, 1 + courant->bd);
    1192                     }
    1193                     else {
    1194                         nouveau = nouvnoeud(PARALLELE, scanne, -1, n);
    1195                         if (DEBUG)
    1196                             printf("Module parallele extrait: %i-%i\n",
    1197                                    1 + debutnoeud->bg, 1 + courant->bd);
    1198                     }
    1199                     extraire(scanne, nouveau, debutfils, nextfils);
    1200                 }
    1201             }
    1202             if (scanne->type == MODULE)
    1203                 scanne->type = PREMIER;
    1204             if (nextfils->suiv == NULL || nextfils->suiv->suiv == NULL
    1205                 || nextfils->suiv->suiv->suiv == NULL)
    1206                 break;
    1207             nextfils = nextfils->suiv;
    1208         }
     1151        nextfils = scanne->fils;        /*on commence au premier fils */
     1152        while (1) {
     1153            courant = nextfils->pointe;
     1154            succourant = nextfils->suiv->pointe;
     1155            if (ps[courant->bd] >= courant->bg
     1156                && ds[courant->bd] <= succourant->bd) {
     1157                /*Des jumeaux!! ->serie ou parallele! */
     1158                /* on va determiner le bloc max de jumeaux consecutifs */
     1159                debutnoeud = courant;
     1160                debutfils = nextfils;
     1161                while (ps[courant->bd] >= courant->bg &&
     1162                       ds[courant->bd] <= succourant->bd &&
     1163                       nextfils->suiv != NULL) {
     1164                    nextfils = nextfils->suiv;
     1165                    courant = nextfils->pointe;
     1166                    if (nextfils->suiv == NULL)
     1167                        break;
     1168                    succourant = nextfils->suiv->pointe;
     1169                }
     1170                /*maintenant on connait la taille du bloc: il va de
     1171                   debutnoeud a courant inclus,
     1172                   et dans la liste des fils de scanne,
     1173                   il s'etend de debutfils a nextfils inclus.
     1174                   On extrait cette liste pour en faire les fils d'un
     1175                   nouveau noeud... sauf si pas la peine! */
     1176                if (debutfils == scanne->fils
     1177                    && nextfils == scanne->lastfils) {
     1178                    /* le noeud scanne etait serie ou parallele */
     1179                  if ( adjii[debutnoeud->bd] !=0)
     1180                        scanne->type = SERIE;
     1181                    else
     1182                        scanne->type = PARALLELE;
     1183                }
     1184                else {
     1185                  if ( adjii[debutnoeud->bd]!=0)
     1186                    /*seule cette ligne distingue G de G' !! */
     1187                    {
     1188                        nouveau = nouvnoeud(SERIE, scanne, -1, n);
     1189                        if (DEBUG)
     1190                            printf("Module serie extrait: %i-%i\n",
     1191                                   1 + debutnoeud->bg, 1 + courant->bd);
     1192                    }
     1193                    else {
     1194                        nouveau = nouvnoeud(PARALLELE, scanne, -1, n);
     1195                        if (DEBUG)
     1196                            printf("Module parallele extrait: %i-%i\n",
     1197                                   1 + debutnoeud->bg, 1 + courant->bd);
     1198                    }
     1199                    extraire(scanne, nouveau, debutfils, nextfils);
     1200                }
     1201            }
     1202            if (scanne->type == MODULE)
     1203                scanne->type = PREMIER;
     1204            if (nextfils->suiv == NULL || nextfils->suiv->suiv == NULL
     1205                || nextfils->suiv->suiv->suiv == NULL)
     1206                break;
     1207            nextfils = nextfils->suiv;
     1208        }
    12091209    }
    12101210    if (DEBUG) {
    1211         printf("Arbre final:\n");
    1212         printnoeud(racine, 0);
     1211        printf("Arbre final:\n");
     1212        printnoeud(racine, 0);
    12131213    }
    12141214    return racine;
    12151215}
     
    12241224      printf("%i : ",i);
    12251225      a=G.G[i];
    12261226      while(a!=NULL)
    1227         {
    1228           printf("%i ", a->s);
    1229           a=a->suiv;
    1230         }
     1227        {
     1228          printf("%i ", a->s);
     1229          a=a->suiv;
     1230        }
    12311231      printf("\n");
    12321232    }
    12331233}
     
    12401240      printf("%i : ",i);
    12411241      a=S[i]->adj;
    12421242      while(a!=NULL)
    1243         {
    1244           printf("%i ", a->pointe->place);
    1245           a=a->suiv;
    1246         }
     1243        {
     1244          printf("%i ", a->pointe->place);
     1245          a=a->suiv;
     1246        }
    12471247      printf("\n");
    12481248    }
    12491249}
     
    12651265/* la fonction principale; qui fait pas grand'chose....*/
    12661266noeud *decomposition_modulaire(graphe G)
    12671267{
    1268  
    1269     sommet **S;                 /* la permutation factorisante */
    1270     noeud *Racine;              /* le futur arbre de decomposition */
    1271  
     1268
     1269    sommet **S;                        /* la permutation factorisante */
     1270    noeud *Racine;                /* le futur arbre de decomposition */
     1271
    12721272    setbuf(stdout,NULL);
    1273    
     1273
    12741274    S = algo1(G);               /* premiere partie: calcul
    1275                                    de la permutation factorisante */
    1276    
     1275                                   de la permutation factorisante */
     1276
    12771277    TrierTous(S,G.n,Calculm(G));/* Trie les listes d'adjacence selon S
    1278                                 */
     1278                                */
    12791279    if(DEBUG)
    12801280      {
    1281         PrintGS(S,G.n);
    1282         PrintS2(S,G.n);
    1283       } 
     1281        PrintGS(S,G.n);
     1282        PrintS2(S,G.n);
     1283      }
    12841284    Racine = algo2(G, S);       /* deuxieme partie: calcul de l'arbre */
    12851285    return Racine;
    12861286}
    1287 
    1288 
    1289 
    1290 
    1291 
    1292 
    1293 
    1294 
    1295 
  • sage/graphs/modular_decomposition/src/random.c

    diff --git a/sage/graphs/modular_decomposition/src/random.c b/sage/graphs/modular_decomposition/src/random.c
    a b  
    2424
    2525#define NIV 20
    2626#define VERBEUX 1
    27  
     27
    2828//extern noeud *decomposition_modulaire(int *,int);
    2929extern void printarbre(noeud *);
    3030/* ppm est la part par million d'arretes voulues dans le graphe */
     
    3939    case PREMIER: C[4*level+2]++; break;
    4040    case FEUILLE: C[4*level+3]++; break;
    4141    }
    42   if(N->type!=FEUILLE) 
     42  if(N->type!=FEUILLE)
    4343    for(F=N->fils;F!=NULL;F=F->suiv)
    4444      compte(F->pointe, level+1, C);
    4545}
     
    6666  for(i=0; i<n; i++)
    6767    {
    6868      if(VERBEUX)
    69         {
    70           printf("Ajacence de %i: ",i+1);
    71           for(a=G.G[i]; a!=NULL; a=a->suiv)
    72             printf("%i ",1+a->s);
    73         }
     69        {
     70          printf("Ajacence de %i: ",i+1);
     71          for(a=G.G[i]; a!=NULL; a=a->suiv)
     72            printf("%i ",1+a->s);
     73        }
    7474      for(j=i+1;j<n;j++)
    75         {
    76         if( (random()%1000000L) < ppm )
    77           {     
    78             // ajoute j a l'adjacence de i
    79             a=(adj *)malloc(sizeof(adj));
    80             a->s=j;
    81             a->suiv=G.G[i];
    82             G.G[i]=a;
    83             // et reciproquement
    84             a=(adj *)malloc(sizeof(adj));
    85             a->s=i;
    86             a->suiv=G.G[j];
    87             G.G[j]=a;
    88             if(VERBEUX)
    89               printf("%i ",j+1);
    90           }
    91         }
     75        {
     76        if( (random()%1000000L) < ppm )
     77          {
     78            // ajoute j a l'adjacence de i
     79            a=(adj *)malloc(sizeof(adj));
     80            a->s=j;
     81            a->suiv=G.G[i];
     82            G.G[i]=a;
     83            // et reciproquement
     84            a=(adj *)malloc(sizeof(adj));
     85            a->s=i;
     86            a->suiv=G.G[j];
     87            G.G[j]=a;
     88            if(VERBEUX)
     89              printf("%i ",j+1);
     90          }
     91        }
    9292      if(VERBEUX)
    93         printf("\n");
     93        printf("\n");
    9494    }
    9595
    9696  // appel de la fonction de decomposition
     
    111111  for(i=1 ; i<NIV ; i++)
    112112    {
    113113      printf("Niveau %i: %i modules (S-P-Pr= %i - %i - %i) et %i feuilles\n",i,
    114        C[4*i]+C[4*i+1]+C[4*i+2], C[4*i], C[4*i+1], C[4*i+2],C[4*i+3]); 
     114       C[4*i]+C[4*i+1]+C[4*i+2], C[4*i], C[4*i+1], C[4*i+2],C[4*i+3]);
    115115      if(i<NIV-1 && C[4*i+4]+C[4*i+5]+C[4*i+6]+C[4*i+7]==0)
    116         break;
     116        break;
    117117    }
    118118  printf("\n");
    119119}
     
    127127  if(narg!=3)
    128128    {
    129129      printf("Decomposition modulaire de graphes \"aleatoires\" \n"
    130              "Donnez en argument:\n"
    131              "le nombre de sommets du graphe\n"
    132              "puis la proportion d'aretes en en millioniemes \n"
    133              "(generateur aleatoire tres primaire)\n"
    134              "Exemple : %s 100 20000\n",arg[0]);
     130             "Donnez en argument:\n"
     131             "le nombre de sommets du graphe\n"
     132             "puis la proportion d'aretes en en millioniemes \n"
     133             "(generateur aleatoire tres primaire)\n"
     134             "Exemple : %s 100 20000\n",arg[0]);
    135135      exit(0);
    136136    }
    137137  n=atoi(arg[1]);
     
    139139  for(i=0;i<4*NIV;i++) C[i]=0;
    140140
    141141  test(n, ppm, C);
    142  
     142
    143143  return 0;
    144144}
    145145
  • sage/libs/linbox/matrix_rational_dense_linbox.cpp

    diff --git a/sage/libs/linbox/matrix_rational_dense_linbox.cpp b/sage/libs/linbox/matrix_rational_dense_linbox.cpp
    a b  
    2121void new_gmp_matrix(DenseMatrix<GMPRationalField>& A, mpq_t** matrix, size_t nrows, size_t ncols) {
    2222    size_t i, j;
    2323    for (i=0; i < nrows; i++) {
    24         for (j=0; j < ncols; j++) {
    25             GMPRationalField::Element t;
    26             t = A.getEntry(i, j);
    27             integer x;
    28             mpz_set(spy.get_mpz(QQ.get_num(x, t)), mpq_numref(matrix[i][j]));
    29             mpz_set(spy.get_mpz(QQ.get_den(x, t)), mpq_denref(matrix[i][j]));
    30             A.setEntry(i, j, t);
    31         }
     24        for (j=0; j < ncols; j++) {
     25            GMPRationalField::Element t;
     26            t = A.getEntry(i, j);
     27            integer x;
     28            mpz_set(spy.get_mpz(QQ.get_num(x, t)), mpq_numref(matrix[i][j]));
     29            mpz_set(spy.get_mpz(QQ.get_den(x, t)), mpq_denref(matrix[i][j]));
     30            A.setEntry(i, j, t);
     31        }
    3232    }
    3333}
    3434 
     
    5353    size_t k;
    5454    k = 19;
    5555    for (size_t i=0; i < nr; i++)
    56         for (size_t j=0; j < nc; j++) {
    57             A.setEntry(i,j,5+i*i+j-j*j+i*i*i+k*k);
    58             k += i*i + j*j + 17;
    59         }
     56        for (size_t j=0; j < nc; j++) {
     57            A.setEntry(i,j,5+i*i+j-j*j+i*i*i+k*k);
     58            k += i*i + j*j + 17;
     59        }
    6060    // new_gmp_matrix(A, matrix, nr, nc);
    6161    cout << "made matrix\n";
    6262    EF.rowEchelon(E, A);
  • sage/libs/mwrank/wrap.cc

    diff --git a/sage/libs/mwrank/wrap.cc b/sage/libs/mwrank/wrap.cc
    a b  
    99
    1010/**************** Miscellaneous functions ****************/
    1111
    12 long mwrank_get_precision() 
     12long mwrank_get_precision()
    1313{
    1414  return decimal_precision();
    1515}
    1616
    17 void mwrank_set_precision(long n) 
     17void mwrank_set_precision(long n)
    1818{
    1919  set_precision(n);
    2020  /*
     
    3030    {n++; set_precision(n);}
    3131}
    3232
    33 void mwrank_initprimes(char* pfilename, int verb) 
     33void mwrank_initprimes(char* pfilename, int verb)
    3434{
    3535  initprimes((const char*)pfilename, verb);
    3636}
     
    6161  return y;
    6262}
    6363
    64 char* bigint_to_str(bigint* x) 
     64char* bigint_to_str(bigint* x)
    6565{
    6666  ostringstream instore;
    6767  instore << (*x);
     
    7171
    7272//////// Curvedata //////////
    7373
    74 struct Curvedata* Curvedata_new(const struct bigint* a1, const struct bigint* a2, 
    75                                 const struct bigint* a3, const struct bigint* a4,
    76                                 const struct bigint* a6, int min_on_init)
     74struct Curvedata* Curvedata_new(const struct bigint* a1, const struct bigint* a2,
     75                                const struct bigint* a3, const struct bigint* a4,
     76                                const struct bigint* a6, int min_on_init)
    7777{
    7878  return new Curvedata(*a1, *a2, *a3, *a4, *a6, min_on_init);
    7979}
     
    9191  return stringstream_to_char(instore);
    9292}
    9393
    94 double Curvedata_silverman_bound(const struct Curvedata* curve) 
     94double Curvedata_silverman_bound(const struct Curvedata* curve)
    9595{
    9696  return silverman_bound(*curve);
    9797}
     
    161161}
    162162
    163163int mw_process(struct Curvedata* curve, struct mw* m,
    164                       const struct bigint* x, const struct bigint* y, 
     164                      const struct bigint* x, const struct bigint* y,
    165165                      const struct bigint* z, int sat)
    166166{
    167167  Point P(*curve, *x, *y, *z);
     
    216216}
    217217
    218218/* Returns index and unsat long array, which user must deallocate */
    219 int mw_saturate(struct mw* m, struct bigint* index, char** unsat, 
     219int mw_saturate(struct mw* m, struct bigint* index, char** unsat,
    220220                       long sat_bd, int odd_primes_only)
    221221{
    222222  vector<long> v;
     
    246246
    247247//////// two_descent //////////
    248248
    249 struct two_descent* two_descent_new(struct Curvedata* curve,  \
    250                                     int verb, int sel,
    251                                     long firstlim, long secondlim,
    252                                     long n_aux, int second_descent)
     249struct two_descent* two_descent_new(struct Curvedata* curve,
     250                                    int verb, int sel,
     251                                    long firstlim, long secondlim,
     252                                    long n_aux, int second_descent)
    253253{
    254254  return new two_descent(curve, verb, sel, firstlim, secondlim, n_aux, second_descent);
    255255}
     
    279279  return p2point_vector_to_str(t->getbasis());
    280280}
    281281
    282 int two_descent_ok(const struct two_descent* t) 
     282int two_descent_ok(const struct two_descent* t)
    283283{
    284284  return t->ok();
    285285}
     
    294294  t->saturate(sat_bd);
    295295}
    296296
    297 char* two_descent_regulator(struct two_descent* t) 
     297char* two_descent_regulator(struct two_descent* t)
    298298{
    299299  bigfloat reg = t->regulator();
    300300  ostringstream instore;
  • sage/libs/mwrank/wrap.h

    diff --git a/sage/libs/mwrank/wrap.h b/sage/libs/mwrank/wrap.h
    a b  
    11#ifdef __cplusplus
    22/* The order here is very important. */
    33#include "eclib/curve.h"
    4 #include "eclib/egr.h"   
     4#include "eclib/egr.h"
    55#include "eclib/descent.h"
    66#include "eclib/points.h"
    77#include "eclib/isogs.h"
    8 #include "eclib/marith.h" 
     8#include "eclib/marith.h"
    99#endif
    1010
    1111#ifdef __cplusplus
    1212#define EXTERN extern "C"
    13 #else 
     13#else
    1414#define EXTERN
    1515#endif
    1616
     
    2727#endif
    2828
    2929#ifdef __cplusplus
    30 extern "C" 
    31 #endif 
     30extern "C"
     31#endif
    3232struct bigint* new_bigint(void);
    3333
    3434EXTERN void del_bigint(struct bigint* x);
     
    4545struct Curvedata;
    4646#endif
    4747
    48 EXTERN struct Curvedata* Curvedata_new(const struct bigint* a1, const struct bigint* a2, 
    49                                        const struct bigint* a3, const struct bigint* a4,
    50                                        const struct bigint* a6, int min_on_init);
     48EXTERN struct Curvedata* Curvedata_new(const struct bigint* a1, const struct bigint* a2,
     49                                       const struct bigint* a3, const struct bigint* a4,
     50                                       const struct bigint* a6, int min_on_init);
    5151
    5252EXTERN void Curvedata_del(struct Curvedata* curve);
    5353
     
    7676EXTERN void mw_del(struct mw* m);
    7777
    7878EXTERN int mw_process(struct Curvedata* curve, struct mw* m,
    79                       const struct bigint* x, const struct bigint* y, 
     79                      const struct bigint* x, const struct bigint* y,
    8080                      const struct bigint* z, int sat);
    8181
    8282EXTERN char* mw_getbasis(struct mw* m);
     
    8686EXTERN int mw_rank(struct mw* m);
    8787
    8888/* Returns index and unsat long array, which user must deallocate */
    89 EXTERN int mw_saturate(struct mw* m, struct bigint* index, char** unsat, 
     89EXTERN int mw_saturate(struct mw* m, struct bigint* index, char** unsat,
    9090                       long sat_bd, int odd_primes_only);
    9191
    9292EXTERN void mw_search(struct mw* m, char* h_lim, int moduli_option, int verb);
     
    9898struct two_descent;
    9999#endif
    100100
    101 EXTERN struct two_descent* two_descent_new(struct Curvedata* curve,  \
    102                                     int verb, int sel,
    103                                     long firstlim, long secondlim,
    104                                     long n_aux, int second_descent);
     101EXTERN struct two_descent* two_descent_new(struct Curvedata* curve,
     102                                    int verb, int sel,
     103                                    long firstlim, long secondlim,
     104                                    long n_aux, int second_descent);
    105105
    106106EXTERN void two_descent_del(struct two_descent* t);
    107107
  • sage/modular/arithgroup/farey.cpp

    diff --git a/sage/modular/arithgroup/farey.cpp b/sage/modular/arithgroup/farey.cpp
    a b  
    143143    vector<long> m;
    144144    for(const_iterator i=gen.begin(); i!=gen.end(); i++) {
    145145      for(const_iterator j=H.begin(); j!=H.end(); j++) {
    146         long q = ((*i)*(*j))%p;
    147         if( find(H.begin(), H.end(), q) == H.end() and
    148             find(m.begin(), m.end(), q) == m.end() ) {
    149           m.push_back(q);
    150         }
     146        long q = ((*i)*(*j))%p;
     147        if( find(H.begin(), H.end(), q) == H.end() and
     148            find(m.begin(), m.end(), q) == m.end() ) {
     149          m.push_back(q);
     150        }
    151151      }
    152152    }
    153153    if( m.size() == 0 ) break;
     
    231231    is_element_general *group = new is_element_general(o);
    232232    // check for user defined SL2Z
    233233    if( group->is_member(SL2Z::S) and
    234         group->is_member(SL2Z::T) ) {
     234        group->is_member(SL2Z::T) ) {
    235235      pairing = vector<int>(2);
    236236      pairing[0] = EVEN;
    237237      pairing[1] = ODD;
     
    247247    }
    248248    // check for index two subgroup
    249249    else if( group->is_member(SL2Z( 0,  1, -1, -1)) and
    250              group->is_member(SL2Z(-1,  1, -1,  0)) ) {
     250             group->is_member(SL2Z(-1,  1, -1,  0)) ) {
    251251      pairing = vector<int>(2);
    252252      pairing[0] = ODD;
    253253      pairing[1] = ODD;
     
    316316    mpq_class largest_diameter(0);
    317317    for(size_t i=0; i<pairing.size(); i++) {
    318318      if( pairing[i] == NO ) {
    319         if( i+1 != pairing.size() ) {
    320           if( i != 0 ) {
    321             mpq_class d = a[i]/b[i] - a[i-1]/b[i-1];
    322             if( d > largest_diameter ) {
    323               largest_diameter = d;
    324               missing_pair = (int)(i);
    325             }
    326           } else {
    327             largest_diameter = infinity;
    328             missing_pair = 0;
    329           }
    330         } else {
    331           largest_diameter = infinity;
    332           missing_pair = (int)(pairing.size()-1);
    333           break;
    334         }
     319        if( i+1 != pairing.size() ) {
     320          if( i != 0 ) {
     321            mpq_class d = a[i]/b[i] - a[i-1]/b[i-1];
     322            if( d > largest_diameter ) {
     323              largest_diameter = d;
     324              missing_pair = (int)(i);
     325            }
     326          } else {
     327            largest_diameter = infinity;
     328            missing_pair = 0;
     329          }
     330        } else {
     331          largest_diameter = infinity;
     332          missing_pair = (int)(pairing.size()-1);
     333          break;
     334        }
    335335      }
    336336    }
    337337    if( missing_pair == -1 ) {
     
    339339    } else {
    340340      mpz_class A, B;
    341341      if( missing_pair+1 == pairing.size() ) {
    342         A = a[missing_pair-1] + 1;
    343         B = b[missing_pair-1] + 0;
     342        A = a[missing_pair-1] + 1;
     343        B = b[missing_pair-1] + 0;
    344344      } else {
    345         if( missing_pair == 0 ) {
    346           A = a[0] - 1;
    347           B = b[0] + 0;
    348         } else {
    349           A = a[missing_pair-1]+a[missing_pair];
    350           B = b[missing_pair-1]+b[missing_pair];
    351         }
     345        if( missing_pair == 0 ) {
     346          A = a[0] - 1;
     347          B = b[0] + 0;
     348        } else {
     349          A = a[missing_pair-1]+a[missing_pair];
     350          B = b[missing_pair-1]+b[missing_pair];
     351        }
    352352      }
    353353      add_term(missing_pair, A/B);
    354354    }
     
    381381  if( pairing[i] == NO ) {
    382382    for(size_t j=0; j<pairing.size(); j++) {
    383383      if( pairing[j] == NO and i != j ) {
    384         vector<int> p(pairing);
    385         p[i] = pairing_max+1;
    386         p[j] = pairing_max+1;
    387         SL2Z C = pairing_matrix(p, i);
    388         if( group->is_member(C) or group->is_member(-C) ) {
    389           pairing_max++;
    390           pairing[i] = pairing_max;
    391           pairing[j] = pairing_max;
    392           return;
    393         }
     384        vector<int> p(pairing);
     385        p[i] = pairing_max+1;
     386        p[j] = pairing_max+1;
     387        SL2Z C = pairing_matrix(p, i);
     388        if( group->is_member(C) or group->is_member(-C) ) {
     389          pairing_max++;
     390          pairing[i] = pairing_max;
     391          pairing[j] = pairing_max;
     392          return;
     393        }
    394394      }
    395395    }
    396396  }
     
    409409    throw(string(__FUNCTION__)+string(": error"));
    410410  } else if( p[i] == EVEN ) {
    411411    return SL2Z(ai1*bi1+ai*bi, -ai*ai-ai1*ai1,
    412                 bi*bi+bi1*bi1, -ai1*bi1-ai*bi);
     412                bi*bi+bi1*bi1, -ai1*bi1-ai*bi);
    413413  } else if( p[i] == ODD ) {
    414414    return SL2Z(ai1*bi1+ai*bi1+ai*bi, -ai*ai-ai*ai1-ai1*ai1,
    415                 bi*bi+bi*bi1+bi1*bi1, -ai1*bi1-ai1*bi-ai*bi);
     415                bi*bi+bi*bi1+bi1*bi1, -ai1*bi1-ai1*bi-ai*bi);
    416416  } else if( p[i] > NO ) {
    417417    const size_t j = paired_side(p, i);
    418418    if( j == 0 ) {
     
    423423      aj = a[j-1]; bj = b[j-1]; aj1 = a[j]; bj1 = b[j];
    424424    }
    425425    return SL2Z(aj1*bi1+aj*bi, -aj*ai-aj1*ai1,
    426                 bj*bi+bj1*bi1, -ai1*bj1-ai*bj);
     426                bj*bi+bj1*bi1, -ai1*bj1-ai*bj);
    427427  }
    428428  return SL2Z::E;
    429429}
     
    518518    size_t i(m), I, J;
    519519    for(;;) {
    520520      if( pairing[i] == NO ) {
    521         I = i;
    522         J = (i==0? pairing.size()-1 : (i-1)%c.size());
     521        I = i;
     522        J = (i==0? pairing.size()-1 : (i-1)%c.size());
    523523      } else {
    524         I = (i+1)%c.size();
    525         J = I;
     524        I = (i+1)%c.size();
     525        J = I;
    526526      }
    527527      if( pairing[I] == ODD or pairing[I] == EVEN ) {
    528         if( c[I] == cusp_class ) {
    529           cusp_class++;
    530           break;
    531         }
    532         c[J] = cusp_class;
    533         i = J;
    534         continue;
     528        if( c[I] == cusp_class ) {
     529          cusp_class++;
     530          break;
     531        }
     532        c[J] = cusp_class;
     533        i = J;
     534        continue;
    535535      } else if( pairing[I] > NO ) {
    536         size_t j;
    537         for(size_t k=0; k<c.size(); k++) {
    538           if( pairing[k] == pairing[I] and k != I ) j = k;
    539         }
    540         if( I != i ) {
    541           if( c[j] == cusp_class ) {
    542             cusp_class++;
    543             break;
    544           }
    545           c[j] = cusp_class;
    546           i = j;
    547           continue;
    548         } else {
    549           if( c[j-1] == cusp_class ) {
    550             cusp_class++;
    551             break;
    552           }
    553           c[j-1] = cusp_class;
    554           i = j-1;
    555           continue;
    556         }
     536        size_t j;
     537        for(size_t k=0; k<c.size(); k++) {
     538          if( pairing[k] == pairing[I] and k != I ) j = k;
     539        }
     540        if( I != i ) {
     541          if( c[j] == cusp_class ) {
     542            cusp_class++;
     543            break;
     544          }
     545          c[j] = cusp_class;
     546          i = j;
     547          continue;
     548        } else {
     549          if( c[j-1] == cusp_class ) {
     550            cusp_class++;
     551            break;
     552          }
     553          c[j-1] = cusp_class;
     554          i = j-1;
     555          continue;
     556        }
    557557      }
    558558    }
    559559  }
     
    566566  for(int i=1; i<=number_of_cusps(); i++) {
    567567    for(size_t j=0; j<cusp_classes.size(); j++) {
    568568      if( cusp_classes[j] == i and cusp_classes[j] != cusp_classes.back() ) {
    569         c.push_back(x[j]);
    570         break;
     569        c.push_back(x[j]);
     570        break;
    571571      }
    572572    }
    573573  }
     
    593593size_t FareySymbol::rank_pi() const {
    594594  if( index() == 2 ) return 1;
    595595  return count_if(pairing.begin(), pairing.end(),
    596                   bind2nd(greater<int>(), 0))/2;
     596                  bind2nd(greater<int>(), 0))/2;
    597597}
    598598
    599599size_t FareySymbol::number_of_cusps() const {
     
    614614    mpq_class cusp_width(0);
    615615    for(size_t j=0; j<cusp_widths.size(); j++) {
    616616      if( cusp_classes[j] == i ) {
    617         cusp_width += cusp_widths[j];
     617        cusp_width += cusp_widths[j];
    618618      }
    619619    }
    620620    width.push_back(cusp_width.get_num());
     
    686686  vector<int> p;
    687687  for(size_t i=0; i<pairing.size(); i++) {
    688688    if( pairing[i] > NO and
    689         p.end() == find(p.begin(), p.end(), pairing[i]) ) {
     689        p.end() == find(p.begin(), p.end(), pairing[i]) ) {
    690690      p.push_back(pairing[i]);
    691691    }
    692692  }
  • sage/modular/arithgroup/farey.hpp

    diff --git a/sage/modular/arithgroup/farey.hpp b/sage/modular/arithgroup/farey.hpp
    a b  
    5353  }
    5454  bool is_member(const SL2Z& V) const {
    5555    return ((V.a()-1) % p == 0 &&
    56             V.c() % p == 0 &&
    57             (V.d()-1) % p == 0);
     56            V.c() % p == 0 &&
     57            (V.d()-1) % p == 0);
    5858  }
    5959};
    6060
     
    6565  }
    6666  bool is_member(const SL2Z& V) const {
    6767    return ((V.a()-1) % p == 0 &&
    68             V.b() % p == 0 &&
    69             V.c() % p == 0 &&
    70             (V.d()-1) % p == 0);
     68            V.b() % p == 0 &&
     69            V.c() % p == 0 &&
     70            (V.d()-1) % p == 0);
    7171  }
    7272};
    7373
  • sage/modular/arithgroup/sl2z.hpp

    diff --git a/sage/modular/arithgroup/sl2z.hpp b/sage/modular/arithgroup/sl2z.hpp
    a b  
    5151
    5252inline
    5353SL2Z::SL2Z(const SL2Z::ElementType& a_, const SL2Z::ElementType& b_,
    54            const SL2Z::ElementType& c_, const SL2Z::ElementType& d_) {
     54           const SL2Z::ElementType& c_, const SL2Z::ElementType& d_) {
    5555  M[0][0] = a_;
    5656  M[0][1] = b_;
    5757  M[1][0] = c_;
     
    8484inline
    8585SL2Z SL2Z::operator*=(const SL2Z& x) {
    8686  return SL2Z(M[0][0]*x.M[0][0] + M[0][1]*x.M[1][0],
    87               M[0][0]*x.M[0][1] + M[0][1]*x.M[1][1],
    88               M[1][0]*x.M[0][0] + M[1][1]*x.M[1][0],
    89               M[1][0]*x.M[0][1] + M[1][1]*x.M[1][1]);
     87              M[0][0]*x.M[0][1] + M[0][1]*x.M[1][1],
     88              M[1][0]*x.M[0][0] + M[1][1]*x.M[1][0],
     89              M[1][0]*x.M[0][1] + M[1][1]*x.M[1][1]);
    9090}
    9191
    9292inline
    9393SL2Z SL2Z::operator/=(const SL2Z& x) {
    9494  return SL2Z( M[0][0]*x.M[1][1] - M[0][1]*x.M[1][0],
    95               -M[0][0]*x.M[0][1] + M[0][1]*x.M[0][0],
    96                M[1][0]*x.M[1][1] - M[1][1]*x.M[1][0],
    97               -M[1][0]*x.M[0][1] + M[1][1]*x.M[0][0]);
     95              -M[0][0]*x.M[0][1] + M[0][1]*x.M[0][0],
     96               M[1][0]*x.M[1][1] - M[1][1]*x.M[1][0],
     97              -M[1][0]*x.M[0][1] + M[1][1]*x.M[0][0]);
    9898}
    9999
    100100inline
     
    110110inline
    111111SL2Z operator*(const SL2Z& x, const SL2Z& y) {
    112112  return SL2Z(x.M[0][0]*y.M[0][0] + x.M[0][1]*y.M[1][0],
    113               x.M[0][0]*y.M[0][1] + x.M[0][1]*y.M[1][1],
    114               x.M[1][0]*y.M[0][0] + x.M[1][1]*y.M[1][0],
    115               x.M[1][0]*y.M[0][1] + x.M[1][1]*y.M[1][1]);
     113              x.M[0][0]*y.M[0][1] + x.M[0][1]*y.M[1][1],
     114              x.M[1][0]*y.M[0][0] + x.M[1][1]*y.M[1][0],
     115              x.M[1][0]*y.M[0][1] + x.M[1][1]*y.M[1][1]);
    116116}
    117117 
    118118inline
    119119SL2Z operator/(const SL2Z& x, const SL2Z& y) {
    120120  return SL2Z( x.M[0][0]*y.M[1][1] - x.M[0][1]*y.M[1][0],
    121               -x.M[0][0]*y.M[0][1] + x.M[0][1]*y.M[0][0],
    122                x.M[1][0]*y.M[1][1] - x.M[1][1]*y.M[1][0],
    123               -x.M[1][0]*y.M[0][1] + x.M[1][1]*y.M[0][0]);
     121              -x.M[0][0]*y.M[0][1] + x.M[0][1]*y.M[0][0],
     122               x.M[1][0]*y.M[1][1] - x.M[1][1]*y.M[1][0],
     123              -x.M[1][0]*y.M[0][1] + x.M[1][1]*y.M[0][0]);
    124124}
    125125 
    126126inline
    127127bool operator==(const SL2Z& x, const SL2Z& y) {
    128128  return (x.M[0][0] == y.M[0][0] and
    129           x.M[0][1] == y.M[0][1] and
    130           x.M[1][0] == y.M[1][0] and
    131           x.M[1][1] == y.M[1][1]);
     129          x.M[0][1] == y.M[0][1] and
     130          x.M[1][0] == y.M[1][0] and
     131          x.M[1][1] == y.M[1][1]);
    132132}
    133133
    134134inline
     
    156156    if( c == ',' ) {
    157157      is >> x.M[0][1] >> c;
    158158      if( c == ';' ) {
    159         is >> x.M[1][0] >> c;
    160         if( c == ',' ) {
    161           is >> x.M[1][1] >> c;
    162           if( c != ']' ) is.clear(std::ios_base::badbit);
    163         } else {
    164           is.clear(std::ios_base::badbit);
    165         }
     159        is >> x.M[1][0] >> c;
     160        if( c == ',' ) {
     161          is >> x.M[1][1] >> c;
     162          if( c != ']' ) is.clear(std::ios_base::badbit);
     163        } else {
     164          is.clear(std::ios_base::badbit);
     165        }
    166166      } else {
    167         is.clear(std::ios_base::badbit);
     167        is.clear(std::ios_base::badbit);
    168168      }
    169169    } else {
    170170      is.clear(std::ios_base::badbit);
  • sage/rings/bernmm/README.txt

    diff --git a/sage/rings/bernmm/README.txt b/sage/rings/bernmm/README.txt
    a b  
    11
    22bernmm: an implementation of the algorithm described in "A multimodular
    3         algorithm for computing Bernoulli numbers", by David Harvey, 2008.
     3        algorithm for computing Bernoulli numbers", by David Harvey, 2008.
    44
    55version 1.1
    66