Compute the sum of powers and the sum of divisors of an integer

AUTHORS

@@ -446,6 +446,8 @@ David Berghaus
 Matthias Gessinger
 
   Graeffe transforms.
+  Sums of powers and sums of divisors of integers.
+  Aliquot sequences.
 
 Erik Postma
 

dev/check_examples.sh

@@ -298,6 +298,23 @@ then
     fi
     echo "PASS"
     exit 0
+elif test "$1" = "aliquot";
+then
+    echo -n "aliquot...."
+    res=$($2/aliquot)
+    if test "$?" != "0";
+    then
+        echo "FAIL"
+        exit 1
+    fi
+    echo "$res" | perl -0ne 'if (/\b117 : 179931895322\b/) { $found=1; last } END { exit !$found }'
+    if test "$?" != "0";
+    then
+        echo "FAIL"
+        exit 2
+    fi
+    echo "PASS"
+    exit 0
 else
     exit 3
 fi

doc/source/fmpz.rst

@@ -959,6 +959,28 @@ Basic arithmetic
 
     Assumes that `m \neq 0`, raises an ``abort`` signal otherwise.
 
+.. function:: void fmpz_sum_powers_horner(fmpz_t f, const fmpz_t g, ulong e)
+              void fmpz_sum_powers_div(fmpz_t f, const fmpz_t g, ulong e)
+	      void fmpz_sum_powers(fmpz_t f, const fmpz_t g, ulong e)
+
+   Sets `f` to the sum of powers of `g`, starting at `g^0 = 1` up to and
+   including `p^e`.
+
+   The `horner` method uses a horner-style method to evaluate the sum
+   iteratively.
+
+   The `div` method uses polynomial division to evaluate the sum with a
+   single exponentiation and a single division.
+
+.. function:: void fmpz_sum_divisors(fmpz_t f, const fmpz_t g)
+
+   Sets `f` to the sum of divisors of `g`, including 1 and `g` itself.
+
+.. function:: void fmpz_sum_divisors_proper(fmpz_t f, const fmpz_t g)
+
+   Sets `f` to the sum of proper divisors of `g`, including 1 but excluding
+   `g` itself.
+
 .. function:: slong fmpz_clog(const fmpz_t x, const fmpz_t b)
               slong fmpz_clog_ui(const fmpz_t x, ulong b)
 

examples/aliquot.c

@@ -0,0 +1,47 @@
+/*
+    Copyright (C) 2024 Matthias Gessinger
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include <flint/fmpz.h>
+#include <flint/profiler.h>
+
+int
+main(int argc, char * argv[])
+{
+    fmpz_t n;
+    ulong num_iterations = 0;
+
+    fmpz_init_set_si(n, 138);
+
+    flint_printf("Computing aliquot sequence\n");
+    flint_printf("%3wu : ", 0);
+    fmpz_print(n);
+    flint_printf("\n");
+
+    TIMEIT_ONCE_START
+
+    while (!fmpz_is_one(n))
+    {
+        fmpz_sum_divisors_proper(n, n);
+        num_iterations++;
+
+        flint_printf("%3wu : ", num_iterations);
+        fmpz_print(n);
+        flint_printf("\n");
+    }
+
+    flint_printf("Sequence terminated after %wu iterations.\n", num_iterations);
+
+    TIMEIT_ONCE_STOP
+    SHOW_MEMORY_USAGE
+
+    fmpz_clear(n);
+    return 0;
+}

src/fmpz.h

@@ -378,6 +378,13 @@ void fmpz_pow_ui(fmpz_t f, const fmpz_t g, ulong exp);
 void fmpz_ui_pow_ui(fmpz_t x, ulong b, ulong e);
 int fmpz_pow_fmpz(fmpz_t a, const fmpz_t b, const fmpz_t e);
 
+void fmpz_sum_powers(fmpz_t f, const fmpz_t g, ulong exp);
+void fmpz_sum_powers_horner(fmpz_t f, const fmpz_t g, ulong exp);
+void fmpz_sum_powers_div(fmpz_t f, const fmpz_t g, ulong exp);
+
+void fmpz_sum_divisors(fmpz_t f, const fmpz_t g);
+void fmpz_sum_divisors_proper(fmpz_t f, const fmpz_t g);
+
 void fmpz_sqrt(fmpz_t f, const fmpz_t g);
 void fmpz_sqrtrem(fmpz_t f, fmpz_t r, const fmpz_t g);
 

src/fmpz/sum_divisors.c

@@ -0,0 +1,75 @@
+/*
+    Copyright (C) 2024 Matthias Gessinger
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "fmpz.h"
+#include "fmpz_factor.h"
+
+void
+fmpz_sum_divisors(fmpz_t f, const fmpz_t g)
+{
+    fmpz_factor_t factors;
+    fmpz * temp;
+    ulong exp;
+    slong i;
+
+    fmpz * p;
+
+    if (fmpz_is_zero(g) || fmpz_is_pm1(g))
+    {
+        fmpz_set(f, g);
+    }
+    else
+    {
+        fmpz_factor_init(factors);
+
+        fmpz_factor(factors, g);
+        i = 0;
+
+        temp = factors->p + i;
+        exp = factors->exp[i];
+
+        fmpz_sum_powers(temp, temp, exp);
+
+        for (i = 1; i < factors->num; i++)
+        {
+            p = factors->p + i;
+            exp = factors->exp[i];
+
+            fmpz_sum_powers(p, p, exp);
+
+            fmpz_mul(temp, temp, p);
+        }
+
+        if (factors->sign == -1)
+        {
+            fmpz_neg(f, temp);
+        }
+        else
+        {
+            fmpz_swap(f, temp);
+        }
+
+        fmpz_factor_clear(factors);
+    }
+}
+
+void
+fmpz_sum_divisors_proper(fmpz_t f, const fmpz_t g)
+{
+    fmpz_t temp;
+    fmpz_init(temp);
+
+    fmpz_sum_divisors(temp, g);
+
+    fmpz_sub(f, temp, g);
+
+    fmpz_clear(temp);
+}

src/fmpz/sum_powers.c

@@ -0,0 +1,92 @@
+/*
+    Copyright (C) 2024 Matthias Gessinger
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "fmpz.h"
+
+void
+fmpz_sum_powers_horner(fmpz_t f, const fmpz_t g, ulong exp)
+{
+    fmpz_t t;
+
+    if (exp == 0 || fmpz_is_zero(g))
+    {
+        fmpz_one(f);
+    }
+    else if (fmpz_is_one(g))
+    {
+        fmpz_set_ui(f, exp + 1);
+    }
+    else
+    {
+        fmpz_init_set_ui(t, 1);
+
+        for (ulong i = 1; i <= exp; i++)
+        {
+            fmpz_mul(t, t, g);
+            fmpz_add_ui(t, t, 1);
+        }
+
+        fmpz_swap(f, t);
+        fmpz_clear(t);
+    }
+}
+
+void
+fmpz_sum_powers_div(fmpz_t f, const fmpz_t g, ulong exp)
+{
+    fmpz_t t;
+
+    if (exp == 0 || fmpz_is_zero(g))
+    {
+        fmpz_one(f);
+    }
+    else if (fmpz_is_one(g))
+    {
+        fmpz_set_ui(f, exp + 1);
+    }
+    else
+    {
+        fmpz_init(t);
+
+        /* t = g - 1 */
+        fmpz_sub_ui(t, g, 1);
+
+        /* f = g^(e + 1) - 1 */
+        fmpz_pow_ui(f, g, exp + 1);
+        fmpz_sub_ui(f, f, 1);
+
+        /* f = (g^(e + 1) - 1) / (g - 1) */
+        fmpz_tdiv_q(f, f, t);
+
+        fmpz_clear(t);
+    }
+}
+
+void
+fmpz_sum_powers(fmpz_t f, const fmpz_t g, ulong exp)
+{
+    if (exp == 0 || fmpz_is_zero(g))
+    {
+        fmpz_one(f);
+    }
+    else if (exp == 1)
+    {
+        fmpz_add_ui(f, g, 1);
+    }
+    else if (exp <= 100)
+    {
+        fmpz_sum_powers_horner(f, g, exp);
+    }
+    else
+    {
+        fmpz_sum_powers_div(f, g, exp);
+    }
+}

src/fmpz/test/main.c

@@ -176,6 +176,10 @@
 #include "t-submul.c"
 #include "t-submul_si.c"
 #include "t-submul_ui.c"
+#include "t-sum_powers_horner.c"
+#include "t-sum_powers_div.c"
+#include "t-sum_powers.c"
+#include "t-sum_divisors.c"
 #include "t-swap.c"
 #include "t-tdiv_q_2exp.c"
 #include "t-tdiv_q.c"
@@ -346,6 +350,10 @@ test_struct tests[] =
     TEST_FUNCTION(fmpz_submul),
     TEST_FUNCTION(fmpz_submul_si),
     TEST_FUNCTION(fmpz_submul_ui),
+    TEST_FUNCTION(fmpz_sum_powers_horner),
+    TEST_FUNCTION(fmpz_sum_powers_div),
+    TEST_FUNCTION(fmpz_sum_powers),
+    TEST_FUNCTION(fmpz_sum_divisors),
     TEST_FUNCTION(fmpz_swap),
     TEST_FUNCTION(fmpz_tdiv_q_2exp),
     TEST_FUNCTION(fmpz_tdiv_q),