Annotation of OpenXM_contrib2/asir2000/engine/up_lm.c, Revision 1.12
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
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 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: *
1.12 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/up_lm.c,v 1.11 2015/08/14 13:51:55 fujimoto Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include <math.h>
52:
53: extern GC_dont_gc;
54: extern struct oEGT eg_chrem,eg_fore,eg_back;
55: extern int debug_up;
56: extern int up_lazy;
1.4 noro 57: extern int current_ff;
1.1 noro 58:
1.6 noro 59: void fft_mulup_lm(UP n1,UP n2,UP *nr)
1.1 noro 60: {
1.12 ! noro 61: ModNum *f1,*f2,*w,*fr;
! 62: ModNum *frarray[1024];
! 63: int modarray[1024];
! 64: int frarray_index = 0;
! 65: N m,m1,m2,lm_mod;
! 66: int d1,d2,dmin,i,mod,root,d,cond,bound;
! 67: UP r;
! 68:
! 69: if ( !n1 || !n2 ) {
! 70: *nr = 0; return;
! 71: }
! 72: d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
! 73: if ( !d1 || !d2 ) {
! 74: mulup(n1,n2,nr); return;
! 75: }
! 76: getmod_lm(&lm_mod);
! 77: if ( !lm_mod )
! 78: error("fft_mulup_lm : current_mod_lm is not set");
! 79: m = ONEN;
! 80: bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;
! 81: f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
! 82: f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
! 83: w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
! 84: for ( i = 0; i < NPrimes; i++ ) {
! 85: FFT_primes(i,&mod,&root,&d);
! 86: if ( (1<<d) < d1+d2+1 )
! 87: continue;
! 88: modarray[frarray_index] = mod;
! 89: frarray[frarray_index++] = fr
! 90: = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
! 91: uptofmarray(mod,n1,f1);
! 92: uptofmarray(mod,n2,f2);
! 93: cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
! 94: if ( cond )
! 95: error("fft_mulup : error in FFT_pol_product");
! 96: STON(mod,m1); muln(m,m1,&m2); m = m2;
! 97: if ( n_bits(m) > bound ) {
! 98: /* GC_dont_gc = 1; */
! 99: crup_lm(frarray,d1+d2,modarray,frarray_index,m,lm_mod,&r);
! 100: uptolmup(r,nr);
! 101: GC_dont_gc = 0;
! 102: return;
! 103: }
! 104: }
! 105: error("fft_mulup : FFT_primes exhausted");
1.1 noro 106: }
107:
1.6 noro 108: void fft_squareup_lm(UP n1,UP *nr)
1.1 noro 109: {
1.12 ! noro 110: ModNum *f1,*w,*fr;
! 111: ModNum *frarray[1024];
! 112: int modarray[1024];
! 113: int frarray_index = 0;
! 114: N m,m1,m2,lm_mod;
! 115: int d1,dmin,i,mod,root,d,cond,bound;
! 116: UP r;
! 117:
! 118: if ( !n1 ) {
! 119: *nr = 0; return;
! 120: }
! 121: d1 = n1->d; dmin = d1;
! 122: if ( !d1 ) {
! 123: mulup(n1,n1,nr); return;
! 124: }
! 125: getmod_lm(&lm_mod);
! 126: if ( !lm_mod )
! 127: error("fft_squareup_lm : current_mod_lm is not set");
! 128: m = ONEN;
! 129: bound = 2*maxblenup(n1)+int_bits(d1)+2;
! 130: f1 = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));
! 131: w = (ModNum *)ALLOCA(6*(1<<int_bits(2*d1+1))*sizeof(ModNum));
! 132: for ( i = 0; i < NPrimes; i++ ) {
! 133: FFT_primes(i,&mod,&root,&d);
! 134: if ( (1<<d) < 2*d1+1 )
! 135: continue;
! 136: modarray[frarray_index] = mod;
! 137: frarray[frarray_index++] = fr
! 138: = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));
! 139: uptofmarray(mod,n1,f1);
! 140: cond = FFT_pol_square(d1,f1,fr,i,w);
! 141: if ( cond )
! 142: error("fft_mulup : error in FFT_pol_product");
! 143: STON(mod,m1); muln(m,m1,&m2); m = m2;
! 144: if ( n_bits(m) > bound ) {
! 145: /* GC_dont_gc = 1; */
! 146: crup_lm(frarray,2*d1,modarray,frarray_index,m,lm_mod,&r);
! 147: uptolmup(r,nr);
! 148: GC_dont_gc = 0;
! 149: return;
! 150: }
! 151: }
! 152: error("fft_squareup : FFT_primes exhausted");
1.1 noro 153: }
154:
1.6 noro 155: void trunc_fft_mulup_lm(UP n1,UP n2,int dbd,UP *nr)
1.1 noro 156: {
1.12 ! noro 157: ModNum *f1,*f2,*fr,*w;
! 158: ModNum *frarray[1024];
! 159: int modarray[1024];
! 160: int frarray_index = 0;
! 161: N m,m1,m2,lm_mod;
! 162: int d1,d2,dmin,i,mod,root,d,cond,bound;
! 163: UP r;
! 164:
! 165: if ( !n1 || !n2 ) {
! 166: *nr = 0; return;
! 167: }
! 168: d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
! 169: if ( !d1 || !d2 ) {
! 170: tmulup(n1,n2,dbd,nr); return;
! 171: }
! 172: getmod_lm(&lm_mod);
! 173: if ( !lm_mod )
! 174: error("trunc_fft_mulup_lm : current_mod_lm is not set");
! 175: m = ONEN;
! 176: bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;
! 177: f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
! 178: f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
! 179: w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
! 180: for ( i = 0; i < NPrimes; i++ ) {
! 181: FFT_primes(i,&mod,&root,&d);
! 182: if ( (1<<d) < d1+d2+1 )
! 183: continue;
! 184:
! 185: modarray[frarray_index] = mod;
! 186: frarray[frarray_index++] = fr
! 187: = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
! 188: uptofmarray(mod,n1,f1);
! 189: uptofmarray(mod,n2,f2);
! 190: cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
! 191: if ( cond )
! 192: error("fft_mulup : error in FFT_pol_product");
! 193: STON(mod,m1); muln(m,m1,&m2); m = m2;
! 194: if ( n_bits(m) > bound ) {
! 195: /* GC_dont_gc = 1; */
! 196: crup_lm(frarray,MIN(dbd-1,d1+d2),modarray,frarray_index,m,lm_mod,&r);
! 197: uptolmup(r,nr);
! 198: GC_dont_gc = 0;
! 199: return;
! 200: }
! 201: }
! 202: error("trunc_fft_mulup : FFT_primes exhausted");
1.1 noro 203: }
204:
1.6 noro 205: void crup_lm(ModNum **f,int d,int *mod,int index,N m,N lm_mod,UP *r)
1.1 noro 206: {
1.12 ! noro 207: double *k;
! 208: double c2;
! 209: unsigned int *w;
! 210: unsigned int zi,au,al;
! 211: UL a;
! 212: int i,j,l,len;
! 213: UP s,s1,u;
! 214: struct oUP c;
! 215: N t,ci,mi,qn;
! 216: unsigned int **sum;
! 217: unsigned int *sum_b;
! 218: Q q;
! 219:
! 220: if ( !lm_mod )
! 221: error("crup_lm : current_mod_lm is not set");
! 222: k = (double *)ALLOCA((d+1)*sizeof(double));
! 223: for ( j = 0; j <= d; j++ )
! 224: k[j] = 0.5;
! 225: up_lazy = 1;
! 226: sum = (unsigned int **)ALLOCA((d+1)*sizeof(unsigned int *));
! 227: len = (int_bits(index)+n_bits(lm_mod)+31)/32+1;
! 228: w = (unsigned int *)ALLOCA((len+1)*sizeof(unsigned int));
! 229: sum_b = (unsigned int *)MALLOC_ATOMIC((d+1)*len*sizeof(unsigned int));
! 230: for ( j = 0; j <= d; j++ ) {
! 231: sum[j] = sum_b+len*j;
! 232: bzero((char *)sum[j],len*sizeof(unsigned int));
! 233: }
! 234: for ( i = 0, s = 0; i < index; i++ ) {
! 235: divin(m,mod[i],&ci);
! 236: zi = (unsigned int)invm((unsigned int)rem(ci,mod[i]),mod[i]);
! 237:
! 238: STON(zi,t); muln(ci,t,&mi);
! 239: divn(mi,lm_mod,&qn,&t);
! 240: if ( t )
! 241: for ( j = 0; j <= d; j++ ) {
! 242: bzero((char *)w,(len+1)*sizeof(unsigned int));
! 243: muln_1(BD(t),PL(t),(unsigned int)f[i][j],w);
! 244: for ( l = PL(t); l >= 0 && !w[l]; l--);
! 245: l++;
! 246: if ( l )
! 247: addarray_to(w,l,sum[j],len);
! 248: }
! 249: c2 = (double)zi/(double)mod[i];
! 250: for ( j = 0; j <= d; j++ )
! 251: k[j] += c2*f[i][j];
! 252: }
! 253: uiarraytoup(sum,len,d,&s);
! 254: GCFREE(sum_b);
1.1 noro 255:
1.12 ! noro 256: u = UPALLOC(d);
! 257: for ( j = 0; j <= d; j++ ) {
1.1 noro 258: #if 1
1.12 ! noro 259: a = (UL)floor(k[j]);
1.11 fujimoto 260: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
1.12 ! noro 261: au = ((unsigned int *)&a)[1];
! 262: al = ((unsigned int *)&a)[0];
1.1 noro 263: #else
1.12 ! noro 264: al = ((unsigned int *)&a)[1];
! 265: au = ((unsigned int *)&a)[0];
1.1 noro 266: #endif
1.12 ! noro 267: if ( au ) {
! 268: NEWQ(q); SGN(q) = 1; NM(q)=NALLOC(2); DN(q)=0;
! 269: PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
! 270: } else
! 271: UTOQ(al,q);
1.1 noro 272: #else
1.12 ! noro 273: al = (int)floor(k[j]); STOQ(al,q);
1.1 noro 274: #endif
1.12 ! noro 275: u->c[j] = (Num)q;
! 276: }
! 277: for ( j = d; j >= 0 && !u->c[j]; j-- );
! 278: if ( j < 0 )
! 279: u = 0;
! 280: else
! 281: u->d = j;
! 282: divn(m,lm_mod,&qn,&t); NTOQ(t,-1,q);
! 283: c.d = 0; c.c[0] = (Num)q;
! 284: mulup(u,&c,&s1);
! 285: up_lazy = 0;
1.1 noro 286:
1.12 ! noro 287: addup(s,s1,r);
1.1 noro 288: }
289:
1.6 noro 290: void fft_rembymulup_special_lm(UP n1,UP n2,UP inv2,UP *nr)
1.1 noro 291: {
1.12 ! noro 292: int d1,d2,d;
! 293: UP r1,t,s,q,u;
1.1 noro 294:
1.12 ! noro 295: if ( !n2 )
! 296: error("rembymulup : division by 0");
! 297: else if ( !n1 || !n2->d )
! 298: *nr = 0;
! 299: else if ( (d1 = n1->d) < (d2 = n2->d) )
! 300: *nr = n1;
! 301: else {
! 302: d = d1-d2;
! 303: reverseup(n1,d1,&t); truncup(t,d+1,&r1);
! 304: trunc_fft_mulup_lm(r1,inv2,d+1,&t);
! 305: truncup(t,d+1,&s);
! 306: reverseup(s,d,&q);
! 307: trunc_fft_mulup_lm(q,n2,d2,&t); truncup(t,d2,&u);
! 308: truncup(n1,d2,&s);
! 309: subup(s,u,nr);
! 310: }
1.1 noro 311: }
312:
1.6 noro 313: void uptolmup(UP n,UP *nr)
1.1 noro 314: {
1.12 ! noro 315: int i,d;
! 316: Q *c;
! 317: LM *cr;
! 318: UP r;
! 319:
! 320: if ( !n )
! 321: *nr = 0;
! 322: else {
! 323: d = n->d; c = (Q *)n->c;
! 324: *nr = r = UPALLOC(d); cr = (LM *)r->c;
! 325: for ( i = 0; i <= d; i++ )
! 326: qtolm(c[i],&cr[i]);
! 327: for ( i = d; i >= 0 && !cr[i]; i-- );
! 328: if ( i < 0 )
! 329: *nr = 0;
! 330: else
! 331: r->d = i;
! 332: }
1.1 noro 333: }
334:
1.6 noro 335: void save_up(UP obj,char *name)
1.1 noro 336: {
1.12 ! noro 337: P p;
! 338: Obj ret;
! 339: NODE n0,n1;
! 340: STRING s;
! 341: void Pbsave();
! 342:
! 343: uptop(obj,&p);
! 344: MKSTR(s,name);
! 345: MKNODE(n1,s,0); MKNODE(n0,p,n1);
! 346: Pbsave(n0,&ret);
1.1 noro 347: }
348:
1.6 noro 349: void hybrid_powermodup(UP f,UP *xp)
1.1 noro 350: {
1.12 ! noro 351: N n;
! 352: UP x,y,t,invf,s;
! 353: int k;
! 354: LM lm;
! 355:
! 356: getmod_lm(&n);
! 357: if ( !n )
! 358: error("hybrid_powermodup : current_mod_lm is not set");
! 359: MKLM(ONEN,lm);
! 360: x = UPALLOC(1); x->d = 1; x->c[1] = (Num)lm;
! 361: y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
! 362:
! 363: reverseup(f,f->d,&t);
! 364: invmodup(t,f->d,&s); uptolmup(s,&invf);
! 365: for ( k = n_bits(n)-1; k >= 0; k-- ) {
! 366: hybrid_squareup(FF_GFP,y,&t);
! 367: hybrid_rembymulup_special(FF_GFP,t,f,invf,&s);
! 368: y = s;
! 369: if ( n->b[k/32] & (1<<(k%32)) ) {
! 370: mulup(y,x,&t);
! 371: remup(t,f,&s);
! 372: y = s;
! 373: }
! 374: }
! 375: *xp = y;
1.1 noro 376: }
377:
1.6 noro 378: void powermodup(UP f,UP *xp)
1.1 noro 379: {
1.12 ! noro 380: N n;
! 381: UP x,y,t,invf,s;
! 382: int k;
! 383: Num c;
! 384:
! 385: if ( !current_ff )
! 386: error("powermodup : current_ff is not set");
! 387: field_order_ff(&n);
! 388: one_ff(&c);
! 389: x = UPALLOC(1); x->d = 1; x->c[1] = c;
! 390: y = UPALLOC(0); y->d = 0; y->c[0] = c;
! 391:
! 392: reverseup(f,f->d,&t);
! 393: invmodup(t,f->d,&s);
! 394: switch ( current_ff ) {
! 395: case FF_GFP:
! 396: case FF_GFPN:
! 397: uptolmup(s,&invf);
! 398: break;
! 399: case FF_GFS:
! 400: case FF_GFSN:
! 401: invf = s; /* XXX */
! 402: break;
! 403: default:
! 404: error("powermodup : not implemented yet");
! 405: }
! 406: for ( k = n_bits(n)-1; k >= 0; k-- ) {
! 407: ksquareup(y,&t);
! 408: rembymulup_special(t,f,invf,&s);
! 409: y = s;
! 410: if ( n->b[k/32] & (1<<(k%32)) ) {
! 411: mulup(y,x,&t);
! 412: remup(t,f,&s);
! 413: y = s;
! 414: }
! 415: }
! 416: *xp = y;
1.1 noro 417: }
418:
419: /* g^d mod f */
420:
1.6 noro 421: void hybrid_generic_powermodup(UP g,UP f,Q d,UP *xp)
1.1 noro 422: {
1.12 ! noro 423: N e;
! 424: UP x,y,t,invf,s;
! 425: int k;
! 426: LM lm;
! 427:
! 428: e = NM(d);
! 429: MKLM(ONEN,lm);
! 430: y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
! 431: remup(g,f,&x);
! 432: if ( !x ) {
! 433: *xp = !d ? y : 0;
! 434: return;
! 435: } else if ( !x->d ) {
! 436: pwrup(x,d,xp);
! 437: return;
! 438: }
! 439: reverseup(f,f->d,&t);
! 440: invmodup(t,f->d,&invf);
! 441: for ( k = n_bits(e)-1; k >= 0; k-- ) {
! 442: hybrid_squareup(FF_GFP,y,&t);
! 443: hybrid_rembymulup_special(FF_GFP,t,f,invf,&s);
! 444: y = s;
! 445: if ( e->b[k/32] & (1<<(k%32)) ) {
! 446: mulup(y,x,&t);
! 447: remup(t,f,&s);
! 448: y = s;
! 449: }
! 450: }
! 451: *xp = y;
1.1 noro 452: }
453:
1.6 noro 454: void generic_powermodup(UP g,UP f,Q d,UP *xp)
1.1 noro 455: {
1.12 ! noro 456: N e;
! 457: UP x,y,t,invf,s;
! 458: int k;
! 459: Num c;
! 460:
! 461: e = NM(d);
! 462: one_ff(&c);
! 463: y = UPALLOC(0); y->d = 0; y->c[0] = c;
! 464: remup(g,f,&x);
! 465: if ( !x ) {
! 466: *xp = !d ? y : 0;
! 467: return;
! 468: } else if ( !x->d ) {
! 469: pwrup(x,d,xp);
! 470: return;
! 471: }
! 472: reverseup(f,f->d,&t);
! 473: invmodup(t,f->d,&invf);
! 474: for ( k = n_bits(e)-1; k >= 0; k-- ) {
! 475: ksquareup(y,&t);
! 476: rembymulup_special(t,f,invf,&s);
! 477: y = s;
! 478: if ( e->b[k/32] & (1<<(k%32)) ) {
! 479: mulup(y,x,&t);
! 480: remup(t,f,&s);
! 481: y = s;
! 482: }
! 483: }
! 484: *xp = y;
1.1 noro 485: }
486:
1.6 noro 487: void hybrid_powertabup(UP f,UP xp,UP *tab)
1.1 noro 488: {
1.12 ! noro 489: UP y,t,invf;
! 490: int i,d;
! 491: LM lm;
! 492:
! 493: d = f->d;
! 494: MKLM(ONEN,lm);
! 495: y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
! 496: tab[0] = y;
! 497: tab[1] = xp;
! 498:
! 499: reverseup(f,f->d,&t);
! 500: invmodup(t,f->d,&invf);
! 501:
! 502: for ( i = 2; i < d; i++ ) {
! 503: if ( debug_up ){
! 504: fprintf(stderr,".");
! 505: }
! 506: hybrid_mulup(FF_GFP,tab[i-1],xp,&t);
! 507: hybrid_rembymulup_special(FF_GFP,t,f,invf,&tab[i]);
! 508: }
1.1 noro 509: }
510:
1.6 noro 511: void powertabup(UP f,UP xp,UP *tab)
1.1 noro 512: {
1.12 ! noro 513: UP y,t,invf;
! 514: int i,d;
! 515: Num c;
! 516:
! 517: d = f->d;
! 518: one_ff(&c);
! 519: y = UPALLOC(0); y->d = 0; y->c[0] = c;
! 520: tab[0] = y;
! 521: tab[1] = xp;
! 522:
! 523: reverseup(f,f->d,&t);
! 524: invmodup(t,f->d,&invf);
! 525:
! 526: for ( i = 2; i < d; i++ ) {
! 527: if ( debug_up ){
! 528: fprintf(stderr,".");
! 529: }
! 530: kmulup(tab[i-1],xp,&t);
! 531: rembymulup_special(t,f,invf,&tab[i]);
! 532: }
1.1 noro 533: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>