Arithmetic derivative

From Rosetta Code
Arithmetic derivative is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The arithmetic derivative of an integer (more specifically, the Lagarias arithmetic derivative) is a function defined for integers, based on prime factorization, by analogy with the product rule for the derivative of a function that is used in mathematical analysis. Accordingly, for natural numbers n, the arithmetic derivative D(n) is defined as follows:

  • D(0) = D(1) = 0.
  • D(p) = 1 for any prime p.
  • D(mn) = D(m)n + mD(n) for any m,n ∈ N. (Leibniz rule for derivatives).

Additionally, for negative integers the arithmetic derivative may be defined as -D(-n) (n < 0).

Examples

D(2) = 1 and D(3) = 1 (both are prime) so if mn = 2 * 3, D(6) = (1)(3) + (1)(2) = 5.

D(9) = D(3)(3) + D(3)(3) = 6

D(27) = D(3)*9 + D(9)*3 = 9 + 18 = 27

D(30) = D(5)(6) + D(6)(5) = 6 + 5 * 5 = 31.

Task

Find and show the arithmetic derivatives for -99 through 100.

Stretch task

Find (the arithmetic derivative of 10^m) then divided by 7, where m is from 1 to 20.

See also


Go

Library: Go-rcu

Using float64 (finessed a little) to avoid the unpleasantness of math/big for the stretch goal. Assumes that int type is 64 bit. <lang go>package main

import (

   "fmt"
   "rcu"

)

func D(n float64) float64 {

   if n < 0 {
       return -D(-n)
   }
   if n < 2 {
       return 0
   }
   var f []int
   if n < 1e19 {
       f = rcu.PrimeFactors(int(n))
   } else {
       g := int(n / 100)
       f = rcu.PrimeFactors(g)
       f = append(f, []int{2, 2, 5, 5}...)
   }
   c := len(f)
   if c == 1 {
       return 1
   }
   if c == 2 {
       return float64(f[0] + f[1])
   }
   d := n / float64(f[0])
   return D(d)*float64(f[0]) + d

}

func main() {

   ad := make([]int, 200)
   for n := -99; n < 101; n++ {
       ad[n+99] = int(D(float64(n)))
   }
   rcu.PrintTable(ad, 10, 4, false)
   fmt.Println()
   pow := 1.0
   for m := 1; m < 21; m++ {
       pow *= 10
       fmt.Printf("D(10^%-2d) / 7 = %.0f\n", m, D(pow)/7)
   }

}</lang>

Output:
 -75  -77   -1 -272  -24  -49  -34  -96  -20 -123 
  -1 -140  -32  -45  -22 -124   -1  -43 -108 -176 
  -1  -71  -18  -80  -55  -39   -1 -156   -1  -59 
 -26  -72   -1  -61  -18 -192  -51  -33   -1  -92 
  -1  -31  -22  -92  -16  -81   -1  -56  -20  -45 
 -14 -112   -1  -25  -39  -48   -1  -41   -1  -68 
 -16  -21   -1  -60  -12  -19  -14  -80   -1  -31 
  -1  -32  -27  -15  -10  -44   -1  -13  -10  -24 
  -1  -21   -1  -32   -8   -9   -1  -16   -1   -7 
  -6  -12   -1   -5   -1   -4   -1   -1    0    0 
   0    1    1    4    1    5    1   12    6    7 
   1   16    1    9    8   32    1   21    1   24 
  10   13    1   44   10   15   27   32    1   31 
   1   80   14   19   12   60    1   21   16   68 
   1   41    1   48   39   25    1  112   14   45 
  20   56    1   81   16   92   22   31    1   92 
   1   33   51  192   18   61    1   72   26   59 
   1  156    1   39   55   80   18   71    1  176 
 108   43    1  124   22   45   32  140    1  123 
  20   96   34   49   24  272    1   77   75  140 

D(10^1 ) / 7 = 1
D(10^2 ) / 7 = 20
D(10^3 ) / 7 = 300
D(10^4 ) / 7 = 4000
D(10^5 ) / 7 = 50000
D(10^6 ) / 7 = 600000
D(10^7 ) / 7 = 7000000
D(10^8 ) / 7 = 80000000
D(10^9 ) / 7 = 900000000
D(10^10) / 7 = 10000000000
D(10^11) / 7 = 110000000000
D(10^12) / 7 = 1200000000000
D(10^13) / 7 = 13000000000000
D(10^14) / 7 = 140000000000000
D(10^15) / 7 = 1500000000000000
D(10^16) / 7 = 16000000000000000
D(10^17) / 7 = 170000000000000000
D(10^18) / 7 = 1800000000000000000
D(10^19) / 7 = 19000000000000000000
D(10^20) / 7 = 200000000000000000000

J

Implementation: <lang J>D=: {{ +/y%q:1>.|y }}"0</lang>

In other words: find the sum of the argument divided by each of the sequence of prime factors of its absolute value (with a special case for zero -- we use the maximum of either 1 or that absolute value when finding the sequence of prime factors).

Task example:

<lang J> D _99+i.20 10 _75 _77 _1 _272 _24 _49 _34 _96 _20 _123

_1 _140 _32  _45 _22 _124  _1  _43 _108 _176
_1  _71 _18  _80 _55  _39  _1 _156   _1  _59

_26 _72 _1 _61 _18 _192 _51 _33 _1 _92

_1  _31 _22  _92 _16  _81  _1  _56  _20  _45

_14 _112 _1 _25 _39 _48 _1 _41 _1 _68 _16 _21 _1 _60 _12 _19 _14 _80 _1 _31

_1  _32 _27  _15 _10  _44  _1  _13  _10  _24
_1  _21  _1  _32  _8   _9  _1  _16   _1   _7
_6  _12  _1   _5  _1   _4  _1   _1    0    0
 0    1   1    4   1    5   1   12    6    7
 1   16   1    9   8   32   1   21    1   24
10   13   1   44  10   15  27   32    1   31
 1   80  14   19  12   60   1   21   16   68
 1   41   1   48  39   25   1  112   14   45
20   56   1   81  16   92  22   31    1   92
 1   33  51  192  18   61   1   72   26   59
 1  156   1   39  55   80  18   71    1  176

108 43 1 124 22 45 32 140 1 123

20   96  34   49  24  272   1   77   75  140</lang>

Also, it seems like it's worth verifying that order of evaluation does not create an ambiguity for the value of D:

<lang J> 15 10 6 + 2 3 5 * D 15 10 6 31 31 31</lang>

Stretch task:

<lang j> (D 10x^1+i.4 5)%7

               1                 20                 300                 4000                 50000
          600000            7000000            80000000            900000000           10000000000
    110000000000      1200000000000      13000000000000      140000000000000      1500000000000000

16000000000000000 170000000000000000 1800000000000000000 19000000000000000000 200000000000000000000</lang>

Julia

<lang ruby>using Primes

D(n) = n < 0 ? -D(-n) : n < 2 ? zero(n) : isprime(n) ? one(n) : typeof(n)(sum(e * n ÷ p for (p, e) in eachfactor(n)))

foreach(p -> print(lpad(p[2], 5), p[1] % 10 == 0 ? "\n" : ""), pairs(map(D, -99:100)))

println() for m in 1:20

   println("D for 10^", rpad(m, 3), "divided by 7 is ", D(Int128(10)^m) ÷ 7)

end

</lang>

Output:
  -75  -77   -1 -272  -24  -49  -34  -96  -20 -123
   -1 -140  -32  -45  -22 -124   -1  -43 -108 -176
   -1  -71  -18  -80  -55  -39   -1 -156   -1  -59
  -26  -72   -1  -61  -18 -192  -51  -33   -1  -92
   -1  -31  -22  -92  -16  -81   -1  -56  -20  -45
  -14 -112   -1  -25  -39  -48   -1  -41   -1  -68
  -16  -21   -1  -60  -12  -19  -14  -80   -1  -31
   -1  -32  -27  -15  -10  -44   -1  -13  -10  -24
   -1  -21   -1  -32   -8   -9   -1  -16   -1   -7
   -6  -12   -1   -5   -1   -4   -1   -1    0    0
    0    1    1    4    1    5    1   12    6    7
    1   16    1    9    8   32    1   21    1   24
   10   13    1   44   10   15   27   32    1   31
    1   80   14   19   12   60    1   21   16   68
    1   41    1   48   39   25    1  112   14   45
   20   56    1   81   16   92   22   31    1   92
    1   33   51  192   18   61    1   72   26   59
    1  156    1   39   55   80   18   71    1  176
  108   43    1  124   22   45   32  140    1  123
   20   96   34   49   24  272    1   77   75  140

D for 10^1  divided by 7 is 1
D for 10^2  divided by 7 is 20
D for 10^3  divided by 7 is 300
D for 10^4  divided by 7 is 4000
D for 10^5  divided by 7 is 50000
D for 10^6  divided by 7 is 600000
D for 10^7  divided by 7 is 7000000
D for 10^8  divided by 7 is 80000000
D for 10^9  divided by 7 is 900000000
D for 10^10 divided by 7 is 10000000000
D for 10^11 divided by 7 is 110000000000
D for 10^12 divided by 7 is 1200000000000
D for 10^13 divided by 7 is 13000000000000
D for 10^14 divided by 7 is 140000000000000
D for 10^15 divided by 7 is 1500000000000000
D for 10^16 divided by 7 is 16000000000000000
D for 10^17 divided by 7 is 170000000000000000
D for 10^18 divided by 7 is 1800000000000000000
D for 10^19 divided by 7 is 19000000000000000000
D for 10^20 divided by 7 is 200000000000000000000

Perl

Translation of: J
Library: ntheory

<lang perl>use v5.36; use bigint; no warnings 'uninitialized'; use List::Util <sum max>; use ntheory 'factor';

sub table ($c, @V) { my $t = $c * (my $w = 2 + length max @V); ( sprintf( ('%'.$w.'d')x@V, @V) ) =~ s/.{1,$t}\K/\n/gr }

sub D ($n) {

   my %f;
   return 0 if $n == 0;
   $f{$_}++ for factor (my $nabs = abs $n);
   my $s = sum map { $nabs * $f{$_} / $_  } keys %f;
   $n > 0 ? $s : -$s;

}

say table 10, map { D $_ } -99 .. 100; say join "\n", map { sprintf('D(10**%-2d) / 7 == ', $_) . D(10**$_) / 7 } 1 .. 20; </lang>

Output:
  -75  -77   -1 -272  -24  -49  -34  -96  -20 -123
   -1 -140  -32  -45  -22 -124   -1  -43 -108 -176
   -1  -71  -18  -80  -55  -39   -1 -156   -1  -59
  -26  -72   -1  -61  -18 -192  -51  -33   -1  -92
   -1  -31  -22  -92  -16  -81   -1  -56  -20  -45
  -14 -112   -1  -25  -39  -48   -1  -41   -1  -68
  -16  -21   -1  -60  -12  -19  -14  -80   -1  -31
   -1  -32  -27  -15  -10  -44   -1  -13  -10  -24
   -1  -21   -1  -32   -8   -9   -1  -16   -1   -7
   -6  -12   -1   -5   -1   -4   -1   -1    0    0
    0    1    1    4    1    5    1   12    6    7
    1   16    1    9    8   32    1   21    1   24
   10   13    1   44   10   15   27   32    1   31
    1   80   14   19   12   60    1   21   16   68
    1   41    1   48   39   25    1  112   14   45
   20   56    1   81   16   92   22   31    1   92
    1   33   51  192   18   61    1   72   26   59
    1  156    1   39   55   80   18   71    1  176
  108   43    1  124   22   45   32  140    1  123
   20   96   34   49   24  272    1   77   75  140

D(10**1 ) / 7 == 1
D(10**2 ) / 7 == 20
D(10**3 ) / 7 == 300
D(10**4 ) / 7 == 4000
D(10**5 ) / 7 == 50000
D(10**6 ) / 7 == 600000
D(10**7 ) / 7 == 7000000
D(10**8 ) / 7 == 80000000
D(10**9 ) / 7 == 900000000
D(10**10) / 7 == 10000000000
D(10**11) / 7 == 110000000000
D(10**12) / 7 == 1200000000000
D(10**13) / 7 == 13000000000000
D(10**14) / 7 == 140000000000000
D(10**15) / 7 == 1500000000000000
D(10**16) / 7 == 16000000000000000
D(10**17) / 7 == 170000000000000000
D(10**18) / 7 == 1800000000000000000
D(10**19) / 7 == 19000000000000000000
D(10**20) / 7 == 200000000000000000000

Phix

with javascript_semantics
include mpfr.e
procedure D(mpz n)
    integer s = mpz_cmp_si(n,0)
    if s<0 then mpz_neg(n,n) end if
    if mpz_cmp_si(n,2)<0 then
        mpz_set_si(n,0)
    else
        sequence f = mpz_prime_factors(n)
        integer c = sum(vslice(f,2)),
                f1 = f[1][1]
        if c=1 then
            mpz_set_si(n,1)
        elsif c=2 then
            mpz_set_si(n,f1 + iff(length(f)=1?f1:f[2][1]))
        else
            assert(mpz_fdiv_q_ui(n,n,f1)=0)
            mpz d = mpz_init_set(n)
            D(n)
            mpz_mul_si(n,n,f1)
            mpz_add(n,n,d)
        end if
        if s<0 then mpz_neg(n,n) end if
    end if
end procedure
 
sequence res = repeat(0,200)
mpz n = mpz_init()
for i=-99 to 100 do
    mpz_set_si(n,i)
    D(n)
    res[i+100] = mpz_get_str(n)
end for
printf(1,"%s\n\n",{join_by(res,1,10," ",fmt:="%4s")})
for m=1 to 20 do
    mpz_ui_pow_ui(n,10,m)
    D(n)
    assert(mpz_fdiv_q_ui(n,n,7)=0)
    printf(1,"D(10^%d)/7 = %s\n",{m,mpz_get_str(n)})
end for
Output:
 -75  -77   -1 -272  -24  -49  -34  -96  -20 -123
  -1 -140  -32  -45  -22 -124   -1  -43 -108 -176
  -1  -71  -18  -80  -55  -39   -1 -156   -1  -59
 -26  -72   -1  -61  -18 -192  -51  -33   -1  -92
  -1  -31  -22  -92  -16  -81   -1  -56  -20  -45
 -14 -112   -1  -25  -39  -48   -1  -41   -1  -68
 -16  -21   -1  -60  -12  -19  -14  -80   -1  -31
  -1  -32  -27  -15  -10  -44   -1  -13  -10  -24
  -1  -21   -1  -32   -8   -9   -1  -16   -1   -7
  -6  -12   -1   -5   -1   -4   -1   -1    0    0
   0    1    1    4    1    5    1   12    6    7
   1   16    1    9    8   32    1   21    1   24
  10   13    1   44   10   15   27   32    1   31
   1   80   14   19   12   60    1   21   16   68
   1   41    1   48   39   25    1  112   14   45
  20   56    1   81   16   92   22   31    1   92
   1   33   51  192   18   61    1   72   26   59
   1  156    1   39   55   80   18   71    1  176
 108   43    1  124   22   45   32  140    1  123
  20   96   34   49   24  272    1   77   75  140

D(10^1)/7 = 1
D(10^2)/7 = 20
D(10^3)/7 = 300
D(10^4)/7 = 4000
D(10^5)/7 = 50000
D(10^6)/7 = 600000
D(10^7)/7 = 7000000
D(10^8)/7 = 80000000
D(10^9)/7 = 900000000
D(10^10)/7 = 10000000000
D(10^11)/7 = 110000000000
D(10^12)/7 = 1200000000000
D(10^13)/7 = 13000000000000
D(10^14)/7 = 140000000000000
D(10^15)/7 = 1500000000000000
D(10^16)/7 = 16000000000000000
D(10^17)/7 = 170000000000000000
D(10^18)/7 = 1800000000000000000
D(10^19)/7 = 19000000000000000000
D(10^20)/7 = 200000000000000000000

Python

<lang python>from sympy.ntheory import factorint

def D(n):

   if n < 0:
       return -D(-n)
   elif n < 2:
       return 0
   else:
       fdict = factorint(n)
       if len(fdict) == 1 and 1 in fdict: # is prime
           return 1
       return sum([n * e // p for p, e in fdict.items()])

for n in range(-99, 101):

   print('{:5}'.format(D(n)), end='\n' if n % 10 == 0 else )

print() for m in range(1, 21):

   print('(D for 10**{}) divided by 7 is {}'.format(m, D(10 ** m) // 7))

</lang>

Output:
  -75  -77   -1 -272  -24  -49  -34  -96  -20 -123
   -1 -140  -32  -45  -22 -124   -1  -43 -108 -176
   -1  -71  -18  -80  -55  -39   -1 -156   -1  -59
  -26  -72   -1  -61  -18 -192  -51  -33   -1  -92
   -1  -31  -22  -92  -16  -81   -1  -56  -20  -45
  -14 -112   -1  -25  -39  -48   -1  -41   -1  -68
  -16  -21   -1  -60  -12  -19  -14  -80   -1  -31
   -1  -32  -27  -15  -10  -44   -1  -13  -10  -24
   -1  -21   -1  -32   -8   -9   -1  -16   -1   -7
   -6  -12   -1   -5   -1   -4   -1   -1    0    0
    0    1    1    4    1    5    1   12    6    7
    1   16    1    9    8   32    1   21    1   24
   10   13    1   44   10   15   27   32    1   31
    1   80   14   19   12   60    1   21   16   68
    1   41    1   48   39   25    1  112   14   45
   20   56    1   81   16   92   22   31    1   92
    1   33   51  192   18   61    1   72   26   59
    1  156    1   39   55   80   18   71    1  176
  108   43    1  124   22   45   32  140    1  123
   20   96   34   49   24  272    1   77   75  140

(D for 10**1) divided by 7 is 1
(D for 10**2) divided by 7 is 20
(D for 10**3) divided by 7 is 300
(D for 10**4) divided by 7 is 4000
(D for 10**5) divided by 7 is 50000
(D for 10**6) divided by 7 is 600000
(D for 10**7) divided by 7 is 7000000
(D for 10**8) divided by 7 is 80000000
(D for 10**9) divided by 7 is 900000000
(D for 10**10) divided by 7 is 10000000000
(D for 10**11) divided by 7 is 110000000000
(D for 10**12) divided by 7 is 1200000000000
(D for 10**13) divided by 7 is 13000000000000
(D for 10**14) divided by 7 is 140000000000000
(D for 10**15) divided by 7 is 1500000000000000
(D for 10**16) divided by 7 is 16000000000000000
(D for 10**17) divided by 7 is 170000000000000000
(D for 10**18) divided by 7 is 1800000000000000000
(D for 10**19) divided by 7 is 19000000000000000000
(D for 10**20) divided by 7 is 200000000000000000000

Raku

<lang perl6>use Prime::Factor;

multi D (0) { 0 } multi D (1) { 0 } multi D ($n where &is-prime) { 1 } multi D ($n where * < 0 ) { -D -$n } multi D ($n) { sum $n.&prime-factors.Bag.map: { $n × .value / .key } }


put (-99 .. 100).map(&D).batch(10)».fmt("%4d").join: "\n";

put ;

put join "\n", (1..20).map: { sprintf "D(10**%-2d) / 7 == %d", $_, D(10**$_) / 7 }</lang>

Output:
 -75  -77   -1 -272  -24  -49  -34  -96  -20 -123
  -1 -140  -32  -45  -22 -124   -1  -43 -108 -176
  -1  -71  -18  -80  -55  -39   -1 -156   -1  -59
 -26  -72   -1  -61  -18 -192  -51  -33   -1  -92
  -1  -31  -22  -92  -16  -81   -1  -56  -20  -45
 -14 -112   -1  -25  -39  -48   -1  -41   -1  -68
 -16  -21   -1  -60  -12  -19  -14  -80   -1  -31
  -1  -32  -27  -15  -10  -44   -1  -13  -10  -24
  -1  -21   -1  -32   -8   -9   -1  -16   -1   -7
  -6  -12   -1   -5   -1   -4   -1   -1    0    0
   0    1    1    4    1    5    1   12    6    7
   1   16    1    9    8   32    1   21    1   24
  10   13    1   44   10   15   27   32    1   31
   1   80   14   19   12   60    1   21   16   68
   1   41    1   48   39   25    1  112   14   45
  20   56    1   81   16   92   22   31    1   92
   1   33   51  192   18   61    1   72   26   59
   1  156    1   39   55   80   18   71    1  176
 108   43    1  124   22   45   32  140    1  123
  20   96   34   49   24  272    1   77   75  140

D(10**1 ) / 7 == 1
D(10**2 ) / 7 == 20
D(10**3 ) / 7 == 300
D(10**4 ) / 7 == 4000
D(10**5 ) / 7 == 50000
D(10**6 ) / 7 == 600000
D(10**7 ) / 7 == 7000000
D(10**8 ) / 7 == 80000000
D(10**9 ) / 7 == 900000000
D(10**10) / 7 == 10000000000
D(10**11) / 7 == 110000000000
D(10**12) / 7 == 1200000000000
D(10**13) / 7 == 13000000000000
D(10**14) / 7 == 140000000000000
D(10**15) / 7 == 1500000000000000
D(10**16) / 7 == 16000000000000000
D(10**17) / 7 == 170000000000000000
D(10**18) / 7 == 1800000000000000000
D(10**19) / 7 == 19000000000000000000
D(10**20) / 7 == 200000000000000000000

Wren

Library: Wren-big
Library: Wren-fmt

As integer arithmetic in Wren is inaccurate above 2^53 we need to use BigInt here. <lang ecmascript>import "./big" for BigInt import "./fmt" for Fmt

var D = Fn.new { |n|

   if (n < 0) return -D.call(-n)
   if (n < 2) return BigInt.zero
   var f = BigInt.primeFactors(n)
   var c = f.count
   if (c == 1) return BigInt.one
   if (c == 2) return f[0] + f[1]
   var d = n / f[0]
   return D.call(d) * f[0] + d

}

var ad = List.filled(200, 0) for (n in -99..100) ad[n+99] = D.call(BigInt.new(n)) Fmt.tprint("$4i", ad, 10) System.print() for (m in 1..20) {

   Fmt.print("D(10^$-2d) / 7 = $i", m, D.call(BigInt.ten.pow(m))/7)

}</lang>

Output:
 -75  -77   -1 -272  -24  -49  -34  -96  -20 -123 
  -1 -140  -32  -45  -22 -124   -1  -43 -108 -176 
  -1  -71  -18  -80  -55  -39   -1 -156   -1  -59 
 -26  -72   -1  -61  -18 -192  -51  -33   -1  -92 
  -1  -31  -22  -92  -16  -81   -1  -56  -20  -45 
 -14 -112   -1  -25  -39  -48   -1  -41   -1  -68 
 -16  -21   -1  -60  -12  -19  -14  -80   -1  -31 
  -1  -32  -27  -15  -10  -44   -1  -13  -10  -24 
  -1  -21   -1  -32   -8   -9   -1  -16   -1   -7 
  -6  -12   -1   -5   -1   -4   -1   -1    0    0 
   0    1    1    4    1    5    1   12    6    7 
   1   16    1    9    8   32    1   21    1   24 
  10   13    1   44   10   15   27   32    1   31 
   1   80   14   19   12   60    1   21   16   68 
   1   41    1   48   39   25    1  112   14   45 
  20   56    1   81   16   92   22   31    1   92 
   1   33   51  192   18   61    1   72   26   59 
   1  156    1   39   55   80   18   71    1  176 
 108   43    1  124   22   45   32  140    1  123 
  20   96   34   49   24  272    1   77   75  140 

D(10^1 ) / 7 = 1
D(10^2 ) / 7 = 20
D(10^3 ) / 7 = 300
D(10^4 ) / 7 = 4000
D(10^5 ) / 7 = 50000
D(10^6 ) / 7 = 600000
D(10^7 ) / 7 = 7000000
D(10^8 ) / 7 = 80000000
D(10^9 ) / 7 = 900000000
D(10^10) / 7 = 10000000000
D(10^11) / 7 = 110000000000
D(10^12) / 7 = 1200000000000
D(10^13) / 7 = 13000000000000
D(10^14) / 7 = 140000000000000
D(10^15) / 7 = 1500000000000000
D(10^16) / 7 = 16000000000000000
D(10^17) / 7 = 170000000000000000
D(10^18) / 7 = 1800000000000000000
D(10^19) / 7 = 19000000000000000000
D(10^20) / 7 = 200000000000000000000