C   88

gcd

Guest on 25th August 2022 01:57:54 PM

  1. /* gnubc program: gcd */
  2. /*
  3.                 NUMBER THEORY ALGORITHMS
  4. 1.  sign(n);
  5. 2.  abs(n);
  6. 3.  mod(a,b), (b > 0): returns a(mod b);
  7. 4.  int(a,b):  returns the integer part of a/b, b nonzero;
  8. 5.  gcd(m,n);
  9. 6.  gcd1(m,n): gcd(m,n)=gcd1(m,n)m+gcd2(m,n)n;
  10. 7.  gcd2(m,n);
  11. 8.  lcm(m,n);
  12. 9.  inv(a,n): returns the inverse of a mod m;
  13. 10. cong(a,b,m): solves ax=b (mod m);
  14. 11. chinese(a,b,m,n): solves x=a(mod m),x=b(mod n);
  15. 12. gcda(m[],n): finds gcd(m[0],...,m[n-1]) and expresses it
  16.     as a linear combination of m[0],...,m[n-1];
  17. 13. lcma(m[],n): finds lcm(m[0],...,m[n-1]);
  18. 14. chinesea(a[],m[],n): solves x=a[i](mod m[i]), i=1,...,n-1;
  19. 15. chineseab(a[],b[],m[],n): solves b[i]*x=a[i](mod m[i]),
  20.     i=1,...,n-1;
  21. 16. mpower(a,b,c): returns a^b (mod c);
  22. 17. exp(a,b): returns a^b;
  23. 18. mthroot(a,b,m): returns the integer part of the mth-root of a/b.
  24.     Here a,b and m are positive integers, m > 1.
  25. 19. mthrootr(a,b,m,r): this gives the mth-root of a/b truncated
  26. 20. min(x,y): this returns x if x<=y, y if x>y.
  27. 21. gcd3(a,b,c): returns gcd(a,b,c)
  28.     to r decimal places, r>=0.
  29. */
  30. /* sign of an integer a */
  31. /* s(a)=1,-1,0, according as a>0,a<0,a=0 */
  32.  
  33. define sign(a){
  34.         if(a>0) return(1)
  35.         if(a<0) return(-1)
  36.         return(0)
  37. }
  38.  
  39. /* absolute value of an integer n */
  40.  
  41. define abs(n){
  42.         if(n>=0) return(n)
  43.         return(-n)
  44. }
  45.  
  46. /*NOTE: in bc we have */
  47. /* a%b=m(a,b) if a>=0 or a<0 and b divides a */
  48. /* a%b=m(a,b)-b if a<0, b>0, a not divisible by b */
  49. /* a/b=[a/b] if a>=0 or a<0 and b divides a */
  50. /* a/b=[a/b]+1 if a<0, b>0, a not divisible by b */
  51. /* a=b(a/b)+a%b */
  52.  
  53. /* mod(a,b)=the least non-negative remainder when an integer
  54.    a is divided by a  positive integer b */
  55.  
  56. define mod(a,b){
  57.         auto c
  58.         c=a%b
  59.         if(a>=0) return(c)
  60.         if(c==0) return(0)
  61.         return(c+b)    
  62. }
  63.  
  64. /* int(a,b)=integer part of a/b, a, b integers, b != 0 */
  65.  
  66. define int(a,b){
  67.         if(b<0){
  68.         a= -a
  69.         b= -b
  70.         }
  71.         return((a-mod(a,b))/b)
  72. }
  73.  
  74. /*  gcd(m,n) for any integers m and n */
  75. /* Euclid's division algorithm is used. */
  76. /* We use gcd(m,n)=gcd(m,|n|) */
  77.  
  78. define gcd(m,n){
  79.         auto a,b,c
  80.         a=abs(m)         /* a=r[0] */
  81.         if(n==0) return(a)
  82.         b=abs(n)         /* b=r[1] */
  83.         c=a%b            /* c=r[2]=r[0] mod(r[1]) */
  84.         while(c>0){
  85.                 a=b
  86.                 b=c
  87.                 c=a%b    /* c=r[j]=r[j-2] mod(r[j-1]) */
  88.         }
  89.         return(b)
  90. }  
  91.  
  92. /* gcd(a,b)=a*h(a,b)+b*k(a,b) for any integers a and b,
  93.  * where h(a,b) and k(a,b) are calculated as follows:
  94.  * Let q[1],...,q[n] be the quotients in Euclid's algorithm.
  95.  * Then the recurrence relation
  96.  * s[0]=1, s[1]=0, s[j]=s[j-2]-q[j-1]*s[j-1] for j=2,...,n
  97.  * is implemented. Then h(a,b)=s[n].
  98.  * The recurrence relation
  99.  * t[0]=0, t[1]=1, t[j]=t[j-2]-q[j-1]*t[j-1] for j=2,...,n
  100.  * is implemented. Then k(a,b)=t[n].
  101.  */
  102.  
  103. /* gcd1(m,n)=h(m,n).
  104.  * gcd2(m,n)=k(m,n).
  105.  */
  106.  
  107. define gcd1(m,n){
  108.         auto a,b,c,h,k,l,q,t
  109.         if(n==0) return(sign(m))
  110.         a=m                /* a=r[0] */
  111.         b=abs(n)           /* b=r[1] */
  112.         c=mod(a,b)         /* c=r[2] */
  113.         if(c==0) return(0)
  114.         l=1                /* l=s[0] */
  115.         k=0                /* k=s[1] */
  116.         while(c>0){
  117.                 q=(a-c)/b  /* q=q[j-1]=[r[j-2]/r[j-1]] */
  118.                 a=b
  119.                 b=c
  120.                 c=mod(a,b) /* c=r[j],r[j-2]=q[j-1]*r[j-1]+r[j]*/
  121.                 h=l-q*k    /* h=s[j],k=s[j-1], l=s[j-2] */
  122.                 l=k
  123.                 k=h
  124.         }        
  125.         return(k)          /* k=s[j] */
  126. }  
  127.  
  128.  
  129. define gcd2(m,n){
  130.         auto a,b,c,h,k,l,q,t,j
  131.         if(n==0) return(0)
  132.         a=m                /* a=r[0] */
  133.         b=abs(n)           /* b=r[1] */
  134.         c=mod(a,b)         /* c=r[2] */
  135.         if(c==0) return(sign(n))
  136.         l=0                /* l=t[0] */
  137.         k=1                /* k=t[1] */
  138. /*      print "     ","   r[0]=",a," t[0]=",0,"\n"
  139.         print "q[1]=",(a-c)/b,"  r[1]=",b," t[1]=",1,"\n" */
  140. /*j=2*/
  141.         while(c>0){
  142.                 q=(a-c)/b  /* q=q[j-1]=[r[j-2]/r[j-1]] */
  143.                 a=b
  144.                 b=c
  145.                 c=mod(a,b) /* c=r[j],r[j-2]=q[j-1]*r[j-1]+r[j]*/
  146.                 h=l-q*k    /* h=t[j],k=t[j-1], l=t[j-2] */
  147. /*              j=j+1
  148.                 print "q[",j-1,"]=",(a-c)/b,"  r[",j-1,"]=",b," t[",j-1,"]=",h,"\n" */
  149.                 l=k
  150.                 k=h
  151.         }
  152.                 q=(a-c)/b
  153. /*              print "        r[",j,"]=",c," t[",j,"]=",l-q*k,"\n" */
  154.         if(n<0) return(-k)
  155.         return(k)          /* k=t[j] */
  156. }  
  157.  
  158. /* inv(a,n) returns the inverse of a (mod n), n>0 if gcd(a,n)=1.
  159.  * uses the fact that |s[j]|<n in general.
  160.  */
  161. define inv(a,n){
  162. auto s
  163. if(gcd(a,n)>1){print "gcd(a,n)>1\n";return(0)}
  164. s=gcd1(a,n)
  165. if(s<0)return(s+n)
  166. return(s)
  167. }
  168. /*  the congruence mx=p(mod n) */
  169.  
  170. define cong(m,p,n){
  171.         auto a,b,t,x,y,z
  172.         a=gcd(m,n)
  173.         if(mod(p,a)>0){
  174.  "the number of solutions is "; return(0)
  175. }
  176.         b=gcd1(m,n)
  177.         y=n/a
  178.         p=p/a
  179.         z=mod(b*p,y)
  180.         print "The solution is x ="
  181.         for(t=0;t<a;t++){print " ",z+t*y,","}
  182.         print " mod ",n,"\n"
  183.         print "or equivalently, x = ",z," mod "; return(y)
  184. }
  185.  
  186. /* the congruence mx=p(mod n)
  187.  * with printout suppressed-used in general form of the Chinese
  188.  * remainder  theorem .
  189.  */
  190. define cong1(m,p,n){
  191.         auto a,b,x,y
  192.         a=gcd(m,n)
  193.         if(mod(p,a)>0){return(-1)}
  194.         b=gcd1(m,n)
  195.         y=n/a
  196.         p=p/a
  197.         x=mod(b*p,y)
  198.         return(x)
  199. }
  200.  
  201. /* the Chinese remainder theorem for the congruences x=a(mod m)
  202.  * and x=b(mod n), m>0, n>0, a and b arbitrary integers.
  203.  * The construction of O. Ore, American Mathematical Monthly,
  204.  * vol.59,pp.365-370,1952, is implemented.
  205.  */
  206.  
  207. define chinese(a,b,m,n) {
  208.         auto c,d,x,y,z
  209.         d = gcd(m,n)
  210.         if(mod(a-b,d)!=0){
  211.                 "the number of solutions is "
  212.                 return(0)
  213.         }
  214.         x= m/d;y=n/d
  215.         z=(m*n)/d
  216.         c = mod(b*gcd1(x,y)*x + a*gcd2(x,y)*y,z)
  217.         print "the solution mod(",z,") is "; return(c)
  218. }
  219.  
  220. /* gcd(m[0],m[1],...,m[n-1])=b[0]*m[0]+...+b[n-1]*m[n-1] */
  221.  
  222. define gcda(m[],n){
  223.         auto a,i,c,d,k,b[],t
  224.         t=n
  225.         n=n-1
  226.         for(i=0;i<=n;i++) b[i]=m[i]
  227.         for(i=1;i<=n;i++) b[i]=gcd(b[i],b[i-1])
  228.         a=b[n]
  229.         b[n]=m[n]
  230.         k=1
  231.         for(i=n;i>=1;i--) {
  232.                 c=b[i];d=b[i-1]
  233.                 b[i]=gcd1(c,d)*k
  234.                 if(i==1) break
  235.                 k=k*gcd2(c,d)
  236.                 b[i-1]=m[i-1]
  237.         }
  238.         b[0]=k*gcd2(c,d)
  239.         for(i=0;i<n;i++){print "(",b[i],")*",m[i],"+"}
  240.         print "(",b[n],")*",m[n]
  241.         print " = "; return(a)
  242. }
  243.  
  244. /* lcm(a,b) for any integers a and b */
  245.  
  246. define lcm(a,b){
  247.         auto g
  248.         g=gcd(a,b)
  249.         if(g==0) return(0)
  250.         return(abs(a*b)/g)
  251. }
  252.  
  253. /* lcm(m[0],m[1],...,m[n-1]) */
  254.  
  255. define lcma(m[],n){
  256.         auto i,b[]
  257.         for(i=0;i<n;i++) b[i]=m[i]
  258.         for(i=1;i<n;i++) b[i]=lcm(b[i],b[i-1])
  259.         return(b[n-1])
  260. }
  261.  
  262. /*
  263.  * a^b (mod c), a,b,c integers, a,b>=0,c>=1
  264.  */
  265.  
  266. define mpower(a,b,c){
  267.         auto x,y,z
  268.         x=a
  269.         y=b
  270.         z=1
  271. /*      print z," ",x," ",y,"\n"*/
  272.         while(y>0){
  273.                 while(y%2==0){
  274.                         y=y/2
  275.                         x=(x*x)%c
  276.                 }
  277.                 y=y-1
  278.                 z=mod(z*x,c)
  279.         }
  280.         return(z)
  281. }
  282.  
  283. /*     
  284.  * the bth power of a, where a is an integer, b a positive integer.
  285.  *
  286.  */
  287. define exp(a,b){
  288.         auto x,y,z
  289.         x=a
  290.         y=b
  291.         z=1
  292.         while(y>0){
  293.                 while(y%2==0){
  294.                         y=y/2
  295.                         x=x*x
  296.                 }
  297.                 y=y-1
  298.                 z=z*x
  299.         }
  300.         return(z)
  301. }
  302.  
  303. /* We use a discrete form of Newton's method to find [x<sup>1/m</sup>],
  304.  * the integer part of the mth root of a positive integer x.
  305.  * (See K.R. Matthews, Computing m-th roots, The College Mathematics Journal 19
  306.  * (1988) 174-176 and H. Lüneburg's book On the Rational Normal Form of
  307.  * Endomorphisms 1987, B.I. Wissenschaftsverlag, Mannheim/Wien/Zürich.)
  308.  * Note that if x is a positive rational, then [x^(1/m)]=[[x]^(1/m)].
  309.  * The answer is given truncated to r decimal places.
  310.  */*/
  311. define mthrootr(a,b,m,r){
  312.         auto s,x,y,i,z,d,f,t
  313.         f=exp(10,r)
  314.         a=a*exp(f,m)
  315.         a=a/b
  316.         if(a==0){
  317.                 if(r==0){
  318.                         "answer = ";return(0)
  319.                 }
  320.                 "answer = 0."
  321.                 for(i=1;i<r;i++)"
  322.                 return(0)
  323.         }
  324.         s=a
  325.         while(s>1){
  326.                 s=s/2
  327.                 i=i+1
  328.         }
  329.         x=exp(2,1+(i/m))
  330.         while(1==1){
  331.                 y=x
  332.                 z=exp(x,m-1)
  333.                 x=x+int(a-x*z,m*z)
  334.                 if(y<=x)break
  335.         }
  336.         if(r==0){
  337.                  " "answer = ";return(y)
  338.         }
  339.         d=y%f
  340.         print "answer = ",y/f, "
  341.         t=r-length(d)
  342.         for(i=1;i<=t;i++){
  343.          print " "
  344.         }
  345.         return(d)
  346. }
  347.  
  348. /* We use a discrete form of Newton's method to find [x<sup>1/m</sup>],
  349. * the integer part of the mth root of a positive integer x.
  350. * (See K.R. Matthews, Computing m-th roots, The College Mathematics Journal 19
  351. * (1988) 174-176 and H. Lüneburg's book On the Rational Normal Form of
  352. * Endomorphisms 1987, B.I. Wissenschaftsverlag, Mannheim/Wien/Zürich.)
  353. * Note that if x is a positive rational, then [x^(1/m)]=[[x]^(1/m)].
  354. */
  355. define mthroot(a,b,m){
  356.         auto s,x,y,i,z,d,j
  357.         a=a/b
  358.         if(a==0)return(0)
  359.         s=a
  360.         while(s>1){
  361.                 s=s/2
  362.                 i=i+1
  363.         }
  364.         x=exp(2,1+(i/m))
  365. /*      print "nt "x[0]="\n""\n"*/
  366.         j=1
  367.         while(1==1){
  368.                 y=x
  369.                 z=exp(x,m-1)
  370.                 x=x+int(a-x*z,m*z/*print "x[",j,"]=",x,"\n"*/n"*/
  371.                 j=j+1
  372.                 if(y<=x) break
  373.         }
  374.         return(y)/* From the bc manual. */./*
  375.  * b(n,m) returns the binomial coefficient, the number of subsets each
  376.  * containing m elements from a set containing n elements.
  377.  */
  378.  */
  379.  
  380. define binomial(n,m){
  381. auto x,j
  382. if(n<m)return(0)
  383. if(m==0)return(1)
  384. x=1
  385. for(j=1;j<=m;j=j+1) x=x*(n-j+1)/j
  386. return(x)/* min(x,y) */) */
  387.  
  388. define min(x,y){
  389.         if(y<x)return(y)
  390.         return(x)
  391. /* the Chinese remainder theorem for n simultaneous congruences
  392.  * x=a[i] (mod m[i]), i=0,1,...,n-1, where m[i]>0 for all i and
  393.  * a[0],...,a[n-1] are arbitrary integers.
  394.  * the construction of O. Ore, American Mathematical Monthly,
  395.  * vol.59,pp.365-370,1952, is implemented.
  396.  */
  397.  /* the Chinese remainder theorem for the congruences x=a(mod m)
  398.  * and x=b(mod n), m>0, n>0, a and b arbitrary integers.
  399.  * The construction of O. Ore, American Mathematical Monthly,
  400.  * vol.59,pp.365-370,1952, is implemented.
  401.  * Same as chinese(), but with printing suppressed, except that -1
  402.  * is returned if here is no solution.
  403.  */
  404.  */
  405.  
  406. define chinese1(a,b,m,n) {
  407.         auto c,d,x,y,z
  408.         d = gcd(m,n)
  409.         if(mod(a-b,d)!=0){
  410.                 return(-1)
  411.         }
  412.         x= m/d;y=n/d
  413.         z=(m*n)/d
  414.         c = mod(b*gcd1(x,y)*x + a*gcd2(x,y)*y,z)
  415.         return(c)
  416. }
  417.  
  418.  
  419. define chinesea(a[],m[],n){
  420.         auto i,m,x
  421.         m=m[0]
  422.         x=a[0]
  423.         for(i=1;i<n;i++){
  424.                 x=chinese1(a[i],x,m[i],m)
  425.                 if(x==-1)break
  426.                 m=lcm(m[i],m)
  427.         }
  428.         if(x==-1){
  429.                 print "there is no solution\n"
  430.         }\
  431.         else{
  432.                 print "the solution mod(") is "is "; return(x)
  433.         }/* the Chinese remainder theorem for n simultaneous congruences
  434.  * a[i]x=b[i] (mod m[i]), i=0,1,...,n-1, where m[i]>0 for all i.
  435.  * a[0],...,a[n-1] and b[0],...,b[n-1] are arbitrary integers.
  436.  */
  437.  */
  438.  
  439. define chineseb(a[],b[],m[],n){
  440.         auto c[],i,m,s,x
  441.         m=m[0]
  442.         x=cong1(a[0],b[0],m[0])
  443.         for(i=1;i<n;i++){
  444.                 s=cong1(a[i],b[i],m[i])
  445.                 if(s==-1){x=-1;break}
  446.                 c[i]=s
  447.                 x=chinese1(c[i],x,m[i],m)
  448.                 if(x==-1)break
  449.                 m=lcm(m[i],m)
  450.         }
  451.         return(x)
  452.  
  453. }
  454.  
  455. define gcd3(a,b,c){
  456. auto t
  457. t=gcd(a,b)
  458. t=gcd(t,c)
  459. return(

Raw Paste


Login or Register to edit or fork this paste. It's free.