![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to lmi.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to lmi.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/lmi.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3: #include "base.h"
! 4: #include "inline.h"
! 5:
! 6: typedef struct oGEN_LM {
! 7: int id;
! 8: N dense;
! 9: N sparse;
! 10: int bit;
! 11: unsigned int lower;
! 12: } *GEN_LM;
! 13:
! 14: void pwrlm0(N,N,N *);
! 15: void gen_simpn(N,N*);
! 16: void gen_simpn_force(N,N*);
! 17:
! 18: int lm_lazy;
! 19:
! 20: static GEN_LM current_mod_lm;
! 21:
! 22: void random_lm(r)
! 23: LM *r;
! 24: {
! 25: N n,t;
! 26:
! 27: if ( !current_mod_lm )
! 28: error("random_lm : current_mod_lm is not set");
! 29: randomn(current_mod_lm->bit+1,&n);
! 30: gen_simpn_force(n,&t);
! 31: MKLM(t,*r);
! 32: }
! 33:
! 34: void ntosparsen(p,bits)
! 35: N p;
! 36: N *bits;
! 37: {
! 38: int l,i,j,nz;
! 39: N r;
! 40: unsigned int *pb,*rb;
! 41: if ( !p )
! 42: *bits = 0;
! 43: else {
! 44: l = n_bits(p);
! 45: pb = BD(p);
! 46: for ( i = l-1, nz = 0; i >= 0; i-- )
! 47: if ( pb[i>>5]&(1<<i&31) ) nz++;
! 48: *bits = r = NALLOC(nz); PL(r) = nz;
! 49: rb = BD(r);
! 50: for ( i = l-1, j = 0; i >= 0; i-- )
! 51: if ( pb[i>>5]&(1<<i&31) )
! 52: rb[j++] = i;
! 53: }
! 54: }
! 55:
! 56: void setmod_lm(p)
! 57: N p;
! 58: {
! 59: int i;
! 60:
! 61: if ( !current_mod_lm )
! 62: current_mod_lm = (GEN_LM)MALLOC(sizeof(struct oGEN_LM));
! 63: bzero((char *)current_mod_lm,sizeof(struct oGEN_LM));
! 64: current_mod_lm->dense = p;
! 65: ntosparsen(p,¤t_mod_lm->sparse);
! 66: current_mod_lm->bit = n_bits(p);
! 67: current_mod_lm->id = UP2_DENSE; /* use ordinary rep. by default */
! 68: if ( !(current_mod_lm->bit%32) ) {
! 69: for ( i = PL(p)-1; i >= 1; i-- )
! 70: if ( BD(p)[i] != 0xffffffff )
! 71: break;
! 72: if ( !i ) {
! 73: current_mod_lm->lower = (unsigned int)((((UL)1)<<32)-((UL)BD(p)[0]));
! 74: current_mod_lm->id = UP2_SPARSE; /* XXX */
! 75: }
! 76: }
! 77: }
! 78:
! 79: void getmod_lm(p)
! 80: N *p;
! 81: {
! 82: if ( !current_mod_lm )
! 83: *p = 0;
! 84: else
! 85: *p = current_mod_lm->dense;
! 86: }
! 87:
! 88: void simplm(n,r)
! 89: LM n;
! 90: LM *r;
! 91: {
! 92: N rem;
! 93:
! 94: if ( !n )
! 95: *r = 0;
! 96: else if ( NID(n) != N_LM )
! 97: *r = n;
! 98: else {
! 99: gen_simpn(n->body,&rem);
! 100: MKLM(rem,*r);
! 101: }
! 102: }
! 103:
! 104: void qtolm(q,l)
! 105: Q q;
! 106: LM *l;
! 107: {
! 108: N rn;
! 109: LM a,b,c;
! 110:
! 111: if ( !q || (OID(q)==O_N && ((NID(q) == N_LM) || (NID(q) == N_GFPN))) ) { /* XXX */
! 112: *l = (LM)q;
! 113: } else if ( OID(q) == O_N && NID(q) == N_Q ) {
! 114: gen_simpn(NM(q),&rn); MKLM(rn,a);
! 115: if ( SGN(q) < 0 ) {
! 116: chsgnlm(a,&c); a = c;
! 117: }
! 118: if ( DN(q) ) {
! 119: gen_simpn_force(DN(q),&rn); MKLM(rn,b);
! 120: if ( !b )
! 121: error("qtolm : invalid denominator");
! 122: divlm(a,b,l);
! 123: } else
! 124: *l = a;
! 125: } else
! 126: error("qtolm : invalid argument");
! 127: }
! 128:
! 129: #define NZLM(a) ((a)&&(NID(a)==N_LM))
! 130:
! 131: void addlm(a,b,c)
! 132: LM a,b;
! 133: LM *c;
! 134: {
! 135: int s;
! 136: N t,t1;
! 137: LM z;
! 138:
! 139: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
! 140: if ( !a )
! 141: *c = b;
! 142: else if ( !b )
! 143: *c = a;
! 144: else {
! 145: addn(a->body,b->body,&t);
! 146: gen_simpn(t,&t1);
! 147: MKLM(t1,*c);
! 148: }
! 149: }
! 150:
! 151: void sublm(a,b,c)
! 152: LM a,b;
! 153: LM *c;
! 154: {
! 155: int s;
! 156: N t,t1,mod;
! 157: LM z;
! 158:
! 159: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
! 160: if ( !b )
! 161: *c = a;
! 162: else if ( !a )
! 163: chsgnlm(b,c);
! 164: else {
! 165: s = subn(a->body,b->body,&t);
! 166: if ( !s )
! 167: *c = 0;
! 168: else {
! 169: MKLM(t,z);
! 170: if ( s < 0 )
! 171: chsgnlm(z,c);
! 172: else
! 173: *c = z;
! 174: }
! 175: }
! 176: }
! 177:
! 178: void mullm(a,b,c)
! 179: LM a,b;
! 180: LM *c;
! 181: {
! 182: N t,q,r;
! 183: LM z;
! 184:
! 185: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
! 186: if ( !a || !b )
! 187: *c = 0;
! 188: else {
! 189: kmuln(a->body,b->body,&t);
! 190: gen_simpn(t,&r);
! 191: MKLM(r,*c);
! 192: }
! 193: }
! 194:
! 195: void divlm(a,b,c)
! 196: LM a,b;
! 197: LM *c;
! 198: {
! 199: LM r;
! 200: Q t,m,i;
! 201: LM z;
! 202:
! 203: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
! 204: if ( !b )
! 205: error("divlm: division by 0");
! 206: else if ( !a )
! 207: *c = 0;
! 208: else {
! 209: NTOQ(b->body,1,t); NTOQ(current_mod_lm->dense,1,m);
! 210: invl(t,m,&i);
! 211: MKLM(NM(i),r);
! 212: mullm(a,r,c);
! 213: }
! 214: }
! 215:
! 216: void chsgnlm(a,c)
! 217: LM a,*c;
! 218: {
! 219: LM t;
! 220: N s,u;
! 221:
! 222: if ( !a )
! 223: *c = a;
! 224: else {
! 225: qtolm((Q)a,&t); a = t;
! 226: gen_simpn_force(a->body,&s);
! 227: subn(current_mod_lm->dense,s,&u);
! 228: MKLM(u,*c);
! 229: }
! 230: }
! 231:
! 232: void pwrlm(a,b,c)
! 233: LM a;
! 234: Q b;
! 235: LM *c;
! 236: {
! 237: LM t;
! 238: N s;
! 239:
! 240: if ( !b ) {
! 241: MKLM(ONEN,*c);
! 242: } else if ( !a )
! 243: *c = 0;
! 244: else {
! 245: qtolm((Q)a,&t); a = t;
! 246: pwrlm0(a->body,NM(b),&s);
! 247: MKLM(s,*c);
! 248: }
! 249: }
! 250:
! 251: void pwrlm0(a,n,c)
! 252: N a,n;
! 253: N *c;
! 254: {
! 255: N n1,t,t1,t2,r1;
! 256: int r;
! 257:
! 258: if ( !n )
! 259: *c = ONEN;
! 260: else if ( UNIN(n) )
! 261: *c = a;
! 262: else {
! 263: r = divin(n,2,&n1);
! 264: pwrlm0(a,n1,&t);
! 265: kmuln(t,t,&t1); gen_simpn(t1,&r1);
! 266: if ( r ) {
! 267: kmuln(r1,a,&t2); gen_simpn(t2,&r1);
! 268: }
! 269: *c = r1;
! 270: }
! 271: }
! 272:
! 273: int cmplm(a,b)
! 274: LM a,b;
! 275: {
! 276: LM z;
! 277:
! 278: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
! 279: if ( !a )
! 280: if ( !b )
! 281: return 0;
! 282: else
! 283: return -1;
! 284: else if ( !b )
! 285: return 1;
! 286: else
! 287: return cmpn(a->body,b->body);
! 288: }
! 289:
! 290: void remn_special(N,N,int,unsigned int ,N *);
! 291:
! 292: void gen_simpn(a,b)
! 293: N a;
! 294: N *b;
! 295: {
! 296: if ( !current_mod_lm )
! 297: error("gen_simpn: current_mod_lm is not set");
! 298: if ( lm_lazy )
! 299: *b = a;
! 300: else if ( current_mod_lm->id == UP2_SPARSE )
! 301: remn_special(a,current_mod_lm->dense,current_mod_lm->bit,current_mod_lm->lower,b);
! 302: else
! 303: remn(a,current_mod_lm->dense,b);
! 304: }
! 305:
! 306: void gen_simpn_force(a,b)
! 307: N a;
! 308: N *b;
! 309: {
! 310: if ( !current_mod_lm )
! 311: error("gen_simpn_force: current_mod_lm is not set");
! 312: else if ( current_mod_lm->id == UP2_SPARSE )
! 313: remn_special(a,current_mod_lm->dense,current_mod_lm->bit,current_mod_lm->lower,b);
! 314: else
! 315: remn(a,current_mod_lm->dense,b);
! 316: }
! 317: