[Help-gsl] greatest common divisor

[Achim Gaedke]

Hi!

I was not able to find a standard UNIX function for the greatest common 
divisor, so I have written one. I put it into the specfunc section, but 
I am not realy sure, if it fits. Also I am not sure, if ldiv is 
available for many other systems. I could replace it by a suitable 
integer division.

Any suggestions are welcome.

Yours, Achim

Index: specfunc/Makefile.am
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/Makefile.am,v
retrieving revision 1.98
diff -u -r1.98 Makefile.am
--- specfunc/Makefile.am        29 Jul 2004 13:11:23 -0000      1.98
+++ specfunc/Makefile.am        5 Jan 2005 01:16:21 -0000
@@ -1,12 +1,12 @@
 noinst_LTLIBRARIES = libgslspecfunc.la 
 
-pkginclude_HEADERS = gsl_sf.h gsl_sf_airy.h gsl_sf_bessel.h gsl_sf_clausen.h 
gsl_sf_coulomb.h gsl_sf_coupling.h gsl_sf_dawson.h gsl_sf_debye.h 
gsl_sf_dilog.h gsl_sf_elementary.h gsl_sf_ellint.h gsl_sf_elljac.h gsl_sf_erf.h 
gsl_sf_exp.h gsl_sf_expint.h gsl_sf_fermi_dirac.h gsl_sf_gamma.h 
gsl_sf_gegenbauer.h gsl_sf_hyperg.h gsl_sf_laguerre.h gsl_sf_lambert.h 
gsl_sf_legendre.h gsl_sf_log.h gsl_sf_pow_int.h gsl_sf_psi.h gsl_sf_result.h 
gsl_sf_synchrotron.h gsl_sf_transport.h gsl_sf_trig.h gsl_sf_zeta.h 
gsl_specfunc.h
+pkginclude_HEADERS = gsl_sf.h gsl_sf_airy.h gsl_sf_bessel.h gsl_sf_clausen.h 
gsl_sf_coulomb.h gsl_sf_coupling.h gsl_sf_dawson.h gsl_sf_debye.h 
gsl_sf_dilog.h gsl_sf_elementary.h gsl_sf_ellint.h gsl_sf_elljac.h gsl_sf_erf.h 
gsl_sf_exp.h gsl_sf_expint.h gsl_sf_fermi_dirac.h gsl_sf_gamma.h gsl_sf_gcd.h 
gsl_sf_gegenbauer.h gsl_sf_hyperg.h gsl_sf_laguerre.h gsl_sf_lambert.h 
gsl_sf_legendre.h gsl_sf_log.h gsl_sf_pow_int.h gsl_sf_psi.h gsl_sf_result.h 
gsl_sf_synchrotron.h gsl_sf_transport.h gsl_sf_trig.h gsl_sf_zeta.h 
gsl_specfunc.h
 
 noinst_HEADERS = bessel_amp_phase.h bessel_olver.h bessel_temme.h bessel.h 
hyperg.h legendre.h eval.h chebyshev.h cheb_eval.c cheb_eval_mode.c check.h 
error.h
 
 INCLUDES= -I$(top_builddir)
 
-libgslspecfunc_la_SOURCES = airy.c airy_der.c airy_zero.c atanint.c bessel.c 
bessel.h bessel_I0.c bessel_I1.c bessel_In.c bessel_Inu.c bessel_J0.c 
bessel_J1.c bessel_Jn.c bessel_Jnu.c bessel_K0.c bessel_K1.c bessel_Kn.c 
bessel_Knu.c bessel_Y0.c bessel_Y1.c bessel_Yn.c bessel_Ynu.c 
bessel_amp_phase.c bessel_amp_phase.h bessel_i.c bessel_j.c bessel_k.c 
bessel_olver.c bessel_temme.c bessel_y.c bessel_zero.c bessel_sequence.c beta.c 
beta_inc.c clausen.c coulomb.c coupling.c coulomb_bound.c dawson.c debye.c 
dilog.c elementary.c ellint.c elljac.c erfc.c exp.c expint.c expint3.c 
fermi_dirac.c gegenbauer.c gamma.c gamma_inc.c hyperg_0F1.c hyperg_2F0.c 
hyperg_1F1.c hyperg_2F1.c hyperg_U.c hyperg.c laguerre.c lambert.c 
legendre_H3d.c legendre_Qn.c legendre_con.c legendre_poly.c log.c poch.c 
pow_int.c psi.c recurse.h result.c shint.c sinint.c synchrotron.c transport.c 
trig.c zeta.c
+libgslspecfunc_la_SOURCES = airy.c airy_der.c airy_zero.c atanint.c bessel.c 
bessel.h bessel_I0.c bessel_I1.c bessel_In.c bessel_Inu.c bessel_J0.c 
bessel_J1.c bessel_Jn.c bessel_Jnu.c bessel_K0.c bessel_K1.c bessel_Kn.c 
bessel_Knu.c bessel_Y0.c bessel_Y1.c bessel_Yn.c bessel_Ynu.c 
bessel_amp_phase.c bessel_amp_phase.h bessel_i.c bessel_j.c bessel_k.c 
bessel_olver.c bessel_temme.c bessel_y.c bessel_zero.c bessel_sequence.c beta.c 
beta_inc.c clausen.c coulomb.c coupling.c coulomb_bound.c dawson.c debye.c 
dilog.c elementary.c ellint.c elljac.c erfc.c exp.c expint.c expint3.c 
fermi_dirac.c gcd.c gegenbauer.c gamma.c gamma_inc.c hyperg_0F1.c hyperg_2F0.c 
hyperg_1F1.c hyperg_2F1.c hyperg_U.c hyperg.c laguerre.c lambert.c 
legendre_H3d.c legendre_Qn.c legendre_con.c legendre_poly.c log.c poch.c 
pow_int.c psi.c recurse.h result.c shint.c sinint.c synchrotron.c transport.c 
trig.c zeta.c
 
 TESTS = $(check_PROGRAMS)
 
Index: specfunc/gsl_sf.h
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/gsl_sf.h,v
retrieving revision 1.24
diff -u -r1.24 gsl_sf.h
--- specfunc/gsl_sf.h   3 Oct 2001 11:37:25 -0000       1.24
+++ specfunc/gsl_sf.h   5 Jan 2005 01:16:22 -0000
@@ -21,6 +21,7 @@
 #include <gsl/gsl_sf_expint.h>
 #include <gsl/gsl_sf_fermi_dirac.h>
 #include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_sf_gcd.h>
 #include <gsl/gsl_sf_gegenbauer.h>
 #include <gsl/gsl_sf_hyperg.h>
 #include <gsl/gsl_sf_laguerre.h>
Index: specfunc/test_sf.c
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/test_sf.c,v
retrieving revision 1.199
diff -u -r1.199 test_sf.c
--- specfunc/test_sf.c  20 Jun 2004 18:41:44 -0000      1.199
+++ specfunc/test_sf.c  5 Jan 2005 01:16:34 -0000
@@ -1716,6 +1716,17 @@
   return s;
 }
 
+
+int test_gcd() {
+  int s=0;
+  if (gcd(1,1)!=1) s+=1;
+  if (gcd(1,1)!=1) s+=1;
+  if (gcd(100,30)!=10) s+=1;
+  if (gcd(30,50)!=10) s+=1;
+  if (gcd(55,57)!=1) s+=1;
+  return s;
+}
+
 int test_results(void)
 {
   int s = 0;
@@ -1772,6 +1783,7 @@
   gsl_test(test_expint(),     "Exponential/Sine/Cosine Integrals");
   gsl_test(test_fermidirac(), "Fermi-Dirac Functions");
   gsl_test(test_gamma(),      "Gamma Functions");
+  gsl_test(test_gcd(),        "GCD Functions");
   gsl_test(test_gegen(),      "Gegenbauer Polynomials");
   gsl_test(test_hyperg(),     "Hypergeometric Functions");
   gsl_test(test_laguerre(),   "Laguerre Polynomials");
Index: specfunc/test_sf.h
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/test_sf.h,v
retrieving revision 1.21
diff -u -r1.21 test_sf.h
--- specfunc/test_sf.h  11 Aug 2003 09:10:07 -0000      1.21
+++ specfunc/test_sf.h  5 Jan 2005 01:16:35 -0000
@@ -92,6 +92,6 @@
 int test_gamma(void);
 int test_hyperg(void);
 int test_legendre(void);
-
+int test_gcd(void);
 
 #endif /* !TEST_SF_H */

/* specfunc/gcd.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  A. Gaedke */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_sf_gcd.h>
#include <gsl/gsl_errno.h>

unsigned long int  gcd(unsigned long int a, unsigned long int b) {
  if (a==0 || b==0) {
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  if (a<b) {
    unsigned long int tmp;
    tmp=a;
    a=b;
    b=tmp;
  }

  while(b!=0) {
    ldiv_t remainder=ldiv(a,b);
    a=b;
    b=remainder.rem;
  }
  return a;
}

/* specfunc/gsl_sf_gcd.h
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  A. Gaedke */

#ifndef __GSL_SF_GCH_H__
#define __GSL_SF_GCH_H__

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/* 
 * calculates the greatest common divisor using the euclidian algorithm
 * 
 * exceptions: GSL_EDOMAIN if a or b is zero
 */


unsigned long int  gcd(unsigned long int a, unsigned long int b);

__END_DECLS

#endif /* __GSL_SF_GCD_H__ */