Changeset e67c6f in git
- Timestamp:
- Sep 30, 2010, 6:23:06 PM (13 years ago)
- Branches:
- (u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a657104b677b4c461d018cbf3204d72d34ad66a9')
- Children:
- efe1b1022b46cc4680baef24f2618a45b809633b
- Parents:
- aa04bce3b7216af7ccfbba57476aa46739bff8e5
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/general.lib
raa04bc re67c6f 15 15 deleteSublist(iv,l); delete entries given by iv from list l 16 16 factorial(n[,../..]); n factorial (=n!) (type int), [type bigint] 17 fibonacci(n [,p]); nth Fibonacci number [char p]17 fibonacci(n); nth Fibonacci number 18 18 kmemory([n[,v]]); active [allocated] memory in kilobyte 19 19 killall(); kill all user-defined variables … … 172 172 "USAGE: binomial(n,k); n,k integers 173 173 RETURN: binomial(n,k); binomial coefficient n choose k 174 @* - of type bigint (computed in characteristic 0)175 NOTE: In any characteristic, binomial(n,k) = coefficient of x^k in (1+x)^n176 SEE ALSO: prime177 174 EXAMPLE: example binomial; shows an example 178 175 " … … 203 200 /////////////////////////////////////////////////////////////////////////////// 204 201 205 static proc binomp (int n, int k, int p)206 //computes binomial coefficient n choose k in basering of char p > 0207 //binomial(n,k) = coefficient of x^k in (1+x)^n.208 //Let n=q*p^j, gcd(q,p)=1, then (1+x)^n = (1 + x^(p^j))^q. We have209 //binomial(n,k)=0 if k!=l*p^j and binomial(n,l*p^j) = binomial(q,l).210 //Do this reduction first. Then, in denominator and numerator211 //of defining formula for binomial coefficient, reduce those factors212 //mod p which are not divisible by p and cancel common factors p. Hence,213 //if n = h*p+r, k=l*p+s, r,s<p, binomial(n,k) = binomial(r,s)*binomial(h,l)214 {215 int l,q,i= 1,n,1;216 number r=1;217 if ( k > n-k )218 { k = n-k;219 }220 if ( k<=0 or k>n) //trivial cases221 { r = (k==0)*r;222 }223 else224 {225 while ( q mod p == 0 )226 { l = l*p;227 q = q div p;228 } //we have now n=q*l, l=p^j, gcd(q,p)=1;229 if (k mod l != 0 )230 { r = 0;231 }232 else233 { l = k div l;234 n = q mod p;235 k = l mod p; //now 0<= k,n <p, use binom0 for n,k236 q = q div p; //recursion for q,l237 l = l div p; //use binomp for q,l238 r = binom0(n,k)*binomp(q,l,p);239 }240 }241 return(r);242 }243 202 /////////////////////////////////////////////////////////////////////////////// 244 203 … … 260 219 example 261 220 { "EXAMPLE:"; echo = 2; 262 factorial(37); ""; //37! (as long integer)221 factorial(37); 263 222 } 264 223 /////////////////////////////////////////////////////////////////////////////// … … 274 233 bigint f,g,h=1,1,1; 275 234 //------------------------------ computation -------------------------------- 276 for (i i=3; ii<=n; ii=ii+1)235 for (int ii=3; ii<=n; ii++) 277 236 { 278 237 h = f+g; f = g; g = h; … … 282 241 example 283 242 { "EXAMPLE:"; echo = 2; 284 fibonacci(42); ""; //f(42) of type string (as long integer) 285 ring r = 2,x,dp; 286 number b = fibonacci(42,2); //f(42) of type number, computed in r 287 b; 243 fibonacci(42); 288 244 } 289 245 ///////////////////////////////////////////////////////////////////////////////
Note: See TracChangeset
for help on using the changeset viewer.