[BACK]Return to lmi.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Annotation of OpenXM_contrib2/asir2000/engine/lmi.c, Revision 1.2

1.2     ! noro        1: /*
        !             2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
        !             3:  * All rights reserved.
        !             4:  *
        !             5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
        !             6:  * non-exclusive and royalty-free license to use, copy, modify and
        !             7:  * redistribute, solely for non-commercial and non-profit purposes, the
        !             8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
        !             9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
        !            10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
        !            11:  * third party developer retains all rights, including but not limited to
        !            12:  * copyrights, in and to the SOFTWARE.
        !            13:  *
        !            14:  * (1) FLL does not grant you a license in any way for commercial
        !            15:  * purposes. You may use the SOFTWARE only for non-commercial and
        !            16:  * non-profit purposes only, such as academic, research and internal
        !            17:  * business use.
        !            18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
        !            19:  * international copyright treaties. If you make copies of the SOFTWARE,
        !            20:  * with or without modification, as permitted hereunder, you shall affix
        !            21:  * to all such copies of the SOFTWARE the above copyright notice.
        !            22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
        !            23:  * shall be made on your publication or presentation in any form of the
        !            24:  * results obtained by use of the SOFTWARE.
        !            25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
        !            26:  * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification
        !            27:  * for such modification or the source code of the modified part of the
        !            28:  * SOFTWARE.
        !            29:  *
        !            30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
        !            31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
        !            32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
        !            33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
        !            34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
        !            35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
        !            36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
        !            37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
        !            38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
        !            39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
        !            40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
        !            41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
        !            42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
        !            43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
        !            44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
        !            45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
        !            46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
        !            47:  *
        !            48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/lmi.c,v 1.1.1.1 1999/12/03 07:39:08 noro Exp $
        !            49: */
1.1       noro       50: #include "ca.h"
                     51: #include "base.h"
                     52: #include "inline.h"
                     53:
                     54: typedef struct oGEN_LM {
                     55:        int id;
                     56:        N dense;
                     57:        N sparse;
                     58:        int bit;
                     59:        unsigned int lower;
                     60: } *GEN_LM;
                     61:
                     62: void pwrlm0(N,N,N *);
                     63: void gen_simpn(N,N*);
                     64: void gen_simpn_force(N,N*);
                     65:
                     66: int lm_lazy;
                     67:
                     68: static GEN_LM current_mod_lm;
                     69:
                     70: void random_lm(r)
                     71: LM *r;
                     72: {
                     73:        N n,t;
                     74:
                     75:        if ( !current_mod_lm )
                     76:                error("random_lm : current_mod_lm is not set");
                     77:        randomn(current_mod_lm->bit+1,&n);
                     78:        gen_simpn_force(n,&t);
                     79:        MKLM(t,*r);
                     80: }
                     81:
                     82: void ntosparsen(p,bits)
                     83: N p;
                     84: N *bits;
                     85: {
                     86:        int l,i,j,nz;
                     87:        N r;
                     88:        unsigned int *pb,*rb;
                     89:        if ( !p )
                     90:                *bits = 0;
                     91:        else {
                     92:                l = n_bits(p);
                     93:                pb = BD(p);
                     94:                for ( i = l-1, nz = 0; i >= 0; i-- )
                     95:                        if ( pb[i>>5]&(1<<i&31) ) nz++;
                     96:                *bits = r = NALLOC(nz); PL(r) = nz;
                     97:                rb = BD(r);
                     98:                for ( i = l-1, j = 0; i >= 0; i-- )
                     99:                        if ( pb[i>>5]&(1<<i&31) )
                    100:                                rb[j++] = i;
                    101:        }
                    102: }
                    103:
                    104: void setmod_lm(p)
                    105: N p;
                    106: {
                    107:        int i;
                    108:
                    109:        if ( !current_mod_lm )
                    110:                current_mod_lm = (GEN_LM)MALLOC(sizeof(struct oGEN_LM));
                    111:        bzero((char *)current_mod_lm,sizeof(struct oGEN_LM));
                    112:        current_mod_lm->dense = p;
                    113:        ntosparsen(p,&current_mod_lm->sparse);
                    114:        current_mod_lm->bit = n_bits(p);
                    115:        current_mod_lm->id = UP2_DENSE; /* use ordinary rep. by default */
                    116:        if ( !(current_mod_lm->bit%32) ) {
                    117:                for ( i = PL(p)-1; i >= 1; i-- )
                    118:                        if ( BD(p)[i] != 0xffffffff )
                    119:                                break;
                    120:                if ( !i ) {
                    121:                        current_mod_lm->lower = (unsigned int)((((UL)1)<<32)-((UL)BD(p)[0]));
                    122:                        current_mod_lm->id = UP2_SPARSE; /* XXX */
                    123:                }
                    124:        }
                    125: }
                    126:
                    127: void getmod_lm(p)
                    128: N *p;
                    129: {
                    130:        if ( !current_mod_lm )
                    131:                *p = 0;
                    132:        else
                    133:                *p = current_mod_lm->dense;
                    134: }
                    135:
                    136: void simplm(n,r)
                    137: LM n;
                    138: LM *r;
                    139: {
                    140:        N rem;
                    141:
                    142:        if ( !n )
                    143:                *r = 0;
                    144:        else if ( NID(n) != N_LM )
                    145:                *r = n;
                    146:        else {
                    147:                gen_simpn(n->body,&rem);
                    148:                MKLM(rem,*r);
                    149:        }
                    150: }
                    151:
                    152: void qtolm(q,l)
                    153: Q q;
                    154: LM *l;
                    155: {
                    156:        N rn;
                    157:        LM a,b,c;
                    158:
                    159:        if ( !q || (OID(q)==O_N && ((NID(q) == N_LM) || (NID(q) == N_GFPN))) ) { /* XXX */
                    160:                *l = (LM)q;
                    161:        } else if ( OID(q) == O_N && NID(q) == N_Q ) {
                    162:                gen_simpn(NM(q),&rn); MKLM(rn,a);
                    163:                if ( SGN(q) < 0 ) {
                    164:                        chsgnlm(a,&c); a = c;
                    165:                }
                    166:                if ( DN(q) ) {
                    167:                        gen_simpn_force(DN(q),&rn); MKLM(rn,b);
                    168:                        if ( !b )
                    169:                                error("qtolm : invalid denominator");
                    170:                        divlm(a,b,l);
                    171:                } else
                    172:                        *l = a;
                    173:        } else
                    174:                error("qtolm : invalid argument");
                    175: }
                    176:
                    177: #define NZLM(a) ((a)&&(NID(a)==N_LM))
                    178:
                    179: void addlm(a,b,c)
                    180: LM a,b;
                    181: LM *c;
                    182: {
                    183:        int s;
                    184:        N t,t1;
                    185:        LM z;
                    186:
                    187:        qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
                    188:        if ( !a )
                    189:                *c = b;
                    190:        else if ( !b )
                    191:                *c = a;
                    192:        else {
                    193:                addn(a->body,b->body,&t);
                    194:                gen_simpn(t,&t1);
                    195:                MKLM(t1,*c);
                    196:        }
                    197: }
                    198:
                    199: void sublm(a,b,c)
                    200: LM a,b;
                    201: LM *c;
                    202: {
                    203:        int s;
                    204:        N t,t1,mod;
                    205:        LM z;
                    206:
                    207:        qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
                    208:        if ( !b )
                    209:                *c = a;
                    210:        else if ( !a )
                    211:                chsgnlm(b,c);
                    212:        else {
                    213:                s = subn(a->body,b->body,&t);
                    214:                if ( !s )
                    215:                        *c = 0;
                    216:                else {
                    217:                        MKLM(t,z);
                    218:                        if ( s < 0 )
                    219:                                chsgnlm(z,c);
                    220:                        else
                    221:                                *c = z;
                    222:                }
                    223:        }
                    224: }
                    225:
                    226: void mullm(a,b,c)
                    227: LM a,b;
                    228: LM *c;
                    229: {
                    230:        N t,q,r;
                    231:        LM z;
                    232:
                    233:        qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
                    234:        if ( !a || !b )
                    235:                *c = 0;
                    236:        else {
                    237:                kmuln(a->body,b->body,&t);
                    238:                gen_simpn(t,&r);
                    239:                MKLM(r,*c);
                    240:        }
                    241: }
                    242:
                    243: void divlm(a,b,c)
                    244: LM a,b;
                    245: LM *c;
                    246: {
                    247:        LM r;
                    248:        Q t,m,i;
                    249:        LM z;
                    250:
                    251:        qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
                    252:        if ( !b )
                    253:                error("divlm: division by 0");
                    254:        else if ( !a )
                    255:                *c = 0;
                    256:        else {
                    257:                NTOQ(b->body,1,t); NTOQ(current_mod_lm->dense,1,m);
                    258:                invl(t,m,&i);
                    259:                MKLM(NM(i),r);
                    260:                mullm(a,r,c);
                    261:        }
                    262: }
                    263:
                    264: void chsgnlm(a,c)
                    265: LM a,*c;
                    266: {
                    267:        LM t;
                    268:        N s,u;
                    269:
                    270:        if ( !a )
                    271:                *c = a;
                    272:        else {
                    273:                qtolm((Q)a,&t); a = t;
                    274:                gen_simpn_force(a->body,&s);
                    275:                subn(current_mod_lm->dense,s,&u);
                    276:                MKLM(u,*c);
                    277:        }
                    278: }
                    279:
                    280: void pwrlm(a,b,c)
                    281: LM a;
                    282: Q b;
                    283: LM *c;
                    284: {
                    285:        LM t;
                    286:        N s;
                    287:
                    288:        if ( !b ) {
                    289:                MKLM(ONEN,*c);
                    290:        } else if ( !a )
                    291:                *c = 0;
                    292:        else {
                    293:                qtolm((Q)a,&t); a = t;
                    294:                pwrlm0(a->body,NM(b),&s);
                    295:                MKLM(s,*c);
                    296:        }
                    297: }
                    298:
                    299: void pwrlm0(a,n,c)
                    300: N a,n;
                    301: N *c;
                    302: {
                    303:        N n1,t,t1,t2,r1;
                    304:        int r;
                    305:
                    306:        if ( !n )
                    307:                *c = ONEN;
                    308:        else if ( UNIN(n) )
                    309:                *c = a;
                    310:        else {
                    311:                r = divin(n,2,&n1);
                    312:                pwrlm0(a,n1,&t);
                    313:                kmuln(t,t,&t1); gen_simpn(t1,&r1);
                    314:                if ( r ) {
                    315:                        kmuln(r1,a,&t2); gen_simpn(t2,&r1);
                    316:                }
                    317:                *c = r1;
                    318:        }
                    319: }
                    320:
                    321: int cmplm(a,b)
                    322: LM a,b;
                    323: {
                    324:        LM z;
                    325:
                    326:        qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
                    327:        if ( !a )
                    328:                if ( !b )
                    329:                        return 0;
                    330:                else
                    331:                        return -1;
                    332:        else if ( !b )
                    333:                        return 1;
                    334:        else
                    335:                return cmpn(a->body,b->body);
                    336: }
                    337:
                    338: void remn_special(N,N,int,unsigned int ,N *);
                    339:
                    340: void gen_simpn(a,b)
                    341: N a;
                    342: N *b;
                    343: {
                    344:        if ( !current_mod_lm )
                    345:                error("gen_simpn: current_mod_lm is not set");
                    346:        if ( lm_lazy )
                    347:                *b = a;
                    348:        else if ( current_mod_lm->id == UP2_SPARSE )
                    349:                remn_special(a,current_mod_lm->dense,current_mod_lm->bit,current_mod_lm->lower,b);
                    350:        else
                    351:                remn(a,current_mod_lm->dense,b);
                    352: }
                    353:
                    354: void gen_simpn_force(a,b)
                    355: N a;
                    356: N *b;
                    357: {
                    358:        if ( !current_mod_lm )
                    359:                error("gen_simpn_force: current_mod_lm is not set");
                    360:        else if ( current_mod_lm->id == UP2_SPARSE )
                    361:                remn_special(a,current_mod_lm->dense,current_mod_lm->bit,current_mod_lm->lower,b);
                    362:        else
                    363:                remn(a,current_mod_lm->dense,b);
                    364: }
                    365:

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>