Programmers: a bit of code golf?

The Game Room
41

You can get some mileage out of the fact that any factorial will have a lot of zeroes at the end, since there are so many powers of ten. You could probably strip those out in the computation of the factorial. That’d still leave you with bignum-sized numbers, though. I’d need to think about whether this could be generalized.

Another trick, of course, is that any large factorial (6! or higher) will be a multiple of 9, so whatever the digit sum is, it’ll also be a multiple of 9. I’m not exactly sure, off the top of my head, how you’d take advantage of that, though.

42

in r5rs:



(define fact 
  (lambda (n) (if (zero? n) 1 (* n (fact (- n 1))))))

(define list-of-digits (lambda (n) (if (> n 10) (cons (modulo n 10) (list-of-digits (truncate (/ n 10)))) (list n))))

(define sum-list (lambda (lst) (apply + lst)))

(define sum-of-fact (lambda (n) (sum-list (list-of-digits ( fact n)))))

and I got 10539 for (sum-of-fact 1000)

43

Yes, the Haskell Integer type is an unbounded integer, whereas Int is an (at least) 32 bit finite precision integer. If you type in an integer literal into the Haskell top level, you’ll get the type 5 :: Num a => a, so literals are interpreted either as an Integer, an Int, or any other type that satisfies the Num constraints depending on context.

44

For what it’s worth, I probably should’ve written the following for slightly better use of tail-recursion (essentially, foldl-ing rather than foldr-ing):



digits 0 = 0
digits n = let (a, b) = divMod n 10 in b : digits a
sumFactorialDigits n = sum $ digits $ product [1..n]

Still, mathematical cleverness is lacking… But perhaps it’s beyond what anyone knows how to do cleverly.

45

Now for extra credit, port these into esoteric languages such as Brainfuck, Befunge, LOLCODE, INTERCAL, or whitespace.

46

There really should be a mathematically clever way of doing that sum-of-digits-of-factorial problem. A crude approximation for the sum of the digits of a number would be S(n) ~ [4.5log[sub]10/sub] (brackets indicating floor), since log[sub]10/sub is the number of digits, and each digit can be anything from 0 to 9. A better estimate would take into account the fact that we know that the last several digits of a factorial are all zero: For n!, the last z digits are all zero, where z = [n/5] + [n/25] + [n/125] + … (there might be a cleaner way to write this, but that’s probably efficient enough for computer code). So the number of digits of a factorial can be approximated by S(n!) ~ [4.5log[sub]10/sub - 4.5z(n)] . We also know that the factorial is a multiple of 9, so we know that the sum of its digits must also be a multiple of 9: This gives us S(n!) ~ 9[0.5log[sub]10/sub - 0.5z(n)] . And finally, of course, we could convert the base 10 log to a natural log, and then apply Sterling’s formula.

This isn’t perfect, of course (for one thing, I really should have taken Benford’s law into account when converting the sum-of-digits to the log), but I’ve a strong hunch that something like this, a little more rigorous, could yield a formula for the answer that would work without need for bignums or the equivalent.

47

How about this problem:

(Stolen from Project Euler)

Consider the fraction, n/d, where n and d are positive integers. If n<d and GreatestCommonFactor(n,d)=1, it is called a reduced proper fraction.

How many elements would be contained in the set of reduced proper fractions for d ≤ 10^9? 10^20?

48

OK, I implemented what I was talking about in the previous post, plus a couple of other refinements. The code (in c) follows (yes, I know I could and should have used ints for some of the variables, but I can never remember the syntax for retyping properly):


#define firstterm -0.81061466795	// 1/2 ln(2pi) - 1
#define ce10 0.434294481904		// convert base e to base 10

double factorialsum(double n)
{
  double lbe, lb10, nfirstdigs, lfirstdigs, firstdigs;
  double sumfirstdigs, curdig, zdigs, fivepow, nmiddigs,result;

  lbe = firstterm + (n+0.5)*log(n+1.0) - n; // Stirling's approximation for ln(n!)
  lb10 = lbe*ce10; 			// log10(n!)
  nfirstdigs = ceil(log10(n) + 1.0);	// Number of reliable digits from Stirling

  for(lfirstdigs = lb10;lfirstdigs > nfirstdigs;lfirstdigs--);
  firstdigs = floor(pow(10,lfirstdigs)); // Finding the first few digits of n!

  sumfirstdigs = 0.0;				// Find the sum of the first few digits
  while(firstdigs){
    curdig = firstdigs - 10.0*floor(firstdigs/10.0);
    sumfirstdigs += curdig;
    firstdigs = (firstdigs-curdig)/10.0;
  }

  zdigs = 0.0;				// Find the number of zero digits
  fivepow = 5.0;
  while(fivepow < n) {
    zdigs += floor(n/fivepow);
    fivepow *= 5.0;
  }
  nmiddigs = floor(lb10) - zdigs - nfirstdigs; // Digits that are neither first nor final zeroes
  result = 9*floor(sumfirstdigs/9.0 + nmiddigs/2.0); // Other digits average 4.5, and the final result must be a multiple of 9
  return(result);
}

This returns 594 for factorialsum(101) (low by 45), and 10431 for factorialsum(1000) (low by 108). Clearly, I’m not there yet. But on the other hand, my method is certainly mathematically sophisticated, gets pretty close, and uses only commonly-available programming tools (i.e., no bignums).

49

What do you mean by “retyping”? (ETA: Oh, nevermind, I guess this means casting…)

I wouldn’t say the problem with the bignum solution is that it do not use only commonly-available programming tools (even if not built-in, any decent, general-purpose programming language will let you easily code bignums up in a few lines [if your programming language does not let you do this, draw the obvious conclusions…]).

The only problems, such as there are, with using bignums to solve the problem of calculating the sum of the digits of n! is that the bignum calculations eventually get expensive, and one might hope for a less expensive solution, and the bignum calculations show no cleverness, and one might hope for a cleverer solution. (Those hopes may or may not be fulfillable). But the “availability” of bignums oughtn’t really be an issue.

50

Here’s the coprime fraction result I came up with. I can’t seem to do better than O(n^3/2), so I havent been able to calculate a result for greater than 10^7.

ETA: number of fractions for 10^7= 30,396,356,427,241



using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Golf
{
    class Straights
    {
        static int[] tmemo;
        public static void Main(string[] args)
        {
            int n = 100000000;
            tmemo = new int[n+1];
            tmemo[1] = 1;
            setprimetotients(n + 1);
           // tmemo[2] = 1;
            long count = 0;
            for (int d = 2; d <= n; d++ )
            {
                count += totient(d);
            }
            Console.Out.WriteLine(count);                    
            Console.In.ReadLine();
        }             
        public static void setprimetotients(int max)
        {
            for (long i = 2; i < max; i++)
            {
                double sq = Math.Sqrt(i);
                int j = 2;
                for (; j <= sq && (i % j !=0); j++) ;
                if (j > sq)
                {
                    long temp = i;
                    while (temp < max)
                    {
                        tmemo[temp] =(int)(temp - (temp / i));
                        temp = temp * i;

                    }
                }
            }
        }
        public static int totient(int x)
        {
            if (tmemo[x] == 0)
            {                            
            double sq = Math.Sqrt(x);
            bool done = false;
            for (int i = 2; i <= sq && !done; i++)
            {
                if (x % i == 0 )
                {                                                         
                    if (gcd(i, x / i) == 1)
                    {
                        tmemo[x] = totient(i) * totient(x / i);
                        done = true;
                    }

                }
            }
            }
            return tmemo[x];
        }
        public static int gcd(int x, int y)
        {           
            int t;
            while (y != 0)
            {
                t = y;
                y = x % y;
                x = t;
            }
            return x;
        }        
    }
    
}
51

I noticed an interesting pattern. I think it can probably be reduced further, but I’m not sure about that.

In Ruby:



module Golf_Frac
    
    # Obviously impractical solution, used for testing small values
    def slow_fc(d) 
        return 0 if d == 1   
        (1..d).inject(0) do |sum, curd|
            (1...curd).each {|n| sum+= 1 if curd.gcd(n) == 1 }
            sum
        end
    end

    # Probably a pattern in here I'm missing  
    def frac_count(d)
        return 0 if d == 1
        return 1 if d == 2
        
        (2...d).inject(d-1) do |sum, n|
            p_factors = uniq_prime_factors(n)
            if p_factors.empty?  # is prime
                sum + d-n - (d-n)/n
            else
                new_val = d-n
                # For each non-prime, start with a count of all numbers between n and d.  Reduce it by one for each prime factor of n that has a multiple in this range.  This overcounts the multiples of two distinct factors, so the count needs to be adjusted upward.  This again undercounts multiples of three factors, so adjust downward.  Continue through all combinations of prime factors, alternately adjusting up and down.  
                for i in (0...(p_factors.size)) do                    
                    sign = i.even? ? -1 : 1  # Note i counts from 0
                    p_factors.choose(i+1) do |factor_set|
                        p_rad = factor_set.inject {|prod,f| prod * f}
                        new_val += sign * ((d-n)/p_rad)
                    end
                end
                sum + new_val
            end
        end
    end

end

Here are the other functions. They could certainly be improved, but I don’t have time for it now.



class Array
    def choose(csize=self.size)
        # Yields each combination of size csize in the Array, as an Array.
        # Does nothing if self is [], or smaller than requested csize.
        #
        # Note that the combinations are produced according to array
        # position, not actual value. If duplicate or permuted results  
        # are not desired, use Array::uniq before calling choose.
        
        raise RangeError, "Choose size must be positive value" if csize <0
        return self if (self.size < csize) || csize == 0       
        self.each_with_index do |e,i| 
            if csize > 1 
                self[i+1..-1].choose(csize-1) do |ra|
                    yield [e] + ra
                end if self[i+1..-1]
            else 
                yield [e]
            end     
        end
        self     
    end   
end
module Golf_Frac

       
    @@primes = []   
    # Unique prime factors of an integer
    
    # Only call this with positive integers > 1.
    # The result is only valid if the @@primes array has been filled, by calling this for all values less than n.  This doesn't do any sorting, either.
    def uniq_prime_factors(n)
        # return [1] if n == 1
        primes_unchecked = @@primes.dup
        upf = []
        while not primes_unchecked.empty?
            p_i = primes_unchecked.pop
            (upf << p_i) if (n % p_i == 0)
        end
        if upf.empty? # n is prime
            @@primes.push(n)
        end
        if upf[0] == n 
            return []
        else
            return upf
        end
    end    
    
    # For debugging
    def show_primes(n)
        @@primes.sort.each do |p|
            break if p > n
            print p
            print ","
        end
        nil
    end 

end
52

One reason I’m reluctant to play this game is that for some of you, C may be the lowest-level language you’ve used extensively. For me it’s the highest-level language. :smack:

For d as large as 10^9 or, especially, 10^20 the Count Proper Fraction problem seems a number-theory problem rather than a programming problem, but since the following relatively brief code gets an exact answer for 10^7 in less than a minute I’ll post it.



typedef unsigned long long Number;

Number findafactor(Number x, Number y)
{

     if (++y < 6) {
          if (x % 2 == 0) return 2;
          if (x % 3 == 0) return 3;
          if (x % 5 == 0) return 5;
          y = 6;
     }
     for (y -= y % 6; y * y < x; y += 6) {
          if (x % (y + 1) == 0) return y + 1;
          if (x % (y + 5) == 0) return y + 5;
     }
     return x;
}

int main()
{
     Number     d, tot, x, y, z, t;

     for (tot = 0, d = 2; d <= 10000000; d++) {
          for (x = z = d, y = 1; y < x; ) {
               y = findafactor(x, y);
               z -= z / y;
               x /= y;
               while ((t = x / y) * y == x)
                    x = t;
          }
          tot += z;
     }
     printf("10^7 --> %Ld
", tot);
     printf("10^9 --> About 10^18 / 3.2898684
");
     printf("10^20 --> About 10^40 / 3.2898684
");
}

I’ll guess the “magic number” 3.2898684 is a simple relative of e, but will leave that as an exercise for the … professor. :cool:

53

I’m guessing you already know the “magic number” is π[sup]2[/sup]/3, since I can’t see where you calculated it from otherwise.

Although you are a bit off in the last digit you provide…

54

No, I just observed that N^2/f(N) was asymptotically constant.

And I guessed e not pi since my intuition has never accepted the mysterious relationships between primes (or factorings) and pi. :smack:

55

Hm, have the last few problems been so hard to excel at that they killed the thread? (I think that describes what happened to me, anyway)

Well, if you like, we can try to dispel the mystery:

The connection here is in two parts: first, realizing that the “magic number” is equal to twice the sum of the -2nd powers of the positive integers (this part does not involve π, just the defining connection between factorings and the multiplicative structure of the positive integers), and then showing that this sum is equal to π[sup]2[/sup]/3 (this part does involve π, obviously, but no longer involves primes/factorings).

The first part is particularly simple: the “magic number” is definitionally the density ratio of arbitrary pairs of positive integers to coprime pairs whose first element is less than their second. Of course, there’s a 50-50 chance of getting the ordering right, so this is just twice the density ratio of arbitrary pairs of positive integers to coprime pairs.

Note that every pair of positive integers is just a scaled-up version of some coprime pair (i.e., every fraction can be put (uniquely) in lowest terms); thus, arbitrary pairs of positive integers break down into “Those which are already in lowest terms”, “Those which are scaled up from lowest terms by scaling both the ‘numerator’ and ‘denominator’ by 2”, “Those scaled up by factors of 3”, “Those scaled up by factors of 4”, etc. Scaling a set of integers by a factor of X causes its density to divide by X; doing this to both the ‘numerator’ and ‘denominator’ in a pair of integers thus divides density by X[sup]2[/sup]. Accordingly, arbitrary pairs of positive integers break down into sets whose density ratio to the coprime pairs are 1/1 [sup]2[/sup], 1/2[sup]2[/sup], 1/3[sup]2[/sup], 1/4[sup]2[/sup], etc. This establishes that our “magic number” is twice the sum 1/1[sup]2[/sup] + 1/2[sup]2[/sup] + 1/3[sup]2[/sup] + 1/4[sup]2[/sup] + …

The first step is now finished. Before moving to the second step, let us take advantage of the symmetry of squaring to note that our magic number can also be viewed as … + 1/(-4)[sup]2[/sup] + 1/(-3)[sup]2[/sup] + 1/(-2)[sup]2[/sup] + 1/(-1)[sup]2[/sup] + 1/1[sup]2[/sup] + 1/2[sup]2[/sup] + 1/3[sup]2[/sup] + 1/4[sup]2[/sup] + …; i.e., the sum of the -2nd powers of the nonzero integers, both positive and negative.

56

Ok, now for that more surprising second part [long story short, the connection to π enters in through the Mercator series, relating the series of reciprocal integers to rates of change of exponential growth, including rotation]:

For any Y, consider the product of Y - 1 with the sum of all the integer powers of Y; expanding this out, each power appears once added and once subtracted, so as to cumulatively cancel out to zero. Accordingly, in some sense*, whenever Y - 1 is invertible (in particular, nonzero), the sum of all the integer powers of Y is zero. [*: For anyone worrying about convergence issues, we can employ, for example, Cesàro summation for our use of this argument]

As a special case of this, letting K be a base of exponential growth and substituting K[sup]x[/sup] in for Y, we have that the sum of K[sup]nx[/sup] over all integers n is zero, over any region where K[sup]x[/sup] - 1 is invertible. Integrating this sum, we find that x + the sum of K[sup]nx[/sup]/(n ln(K)) over all nonzero integers n is constant over such a region; equivalently, ln(K) * x + the sum of K[sup]nx[/sup]/n over nonzero integers n is constant. [This is the key fact we will exploit to calculate our magic number; note the presence of the reciprocal integers in the coefficients of the sum]

To calculate this constant value, we would like to plug in a convenient value for x; unfortunately, the most convenient choice (choosing x such that K[sup]x[/sup] = 1, under which the terms of our infinite series for opposite nonzero values of n cancel out, leaving only ln(K) times our choice of x) is not available to us, as this would make K[sup]x[/sup] - 1 zero. But we can achieve pretty much the same cancellation effect with another choice: if K[sup]c[/sup] + 1 = 0, then our sum cancels out to simply ln(K) * c in the same way.

All of this was phrased in terms of arbitrary K so far, but in particular, making the convenient choice of K such that K[sup]x[/sup] is the operation of rotation by x many half revolutions, we have that K[sup]1[/sup] + 1 = 0, so that ln(K) * x + the sum of K[sup]nx[/sup]/n over all nonzero integers n is constantly ln(K) * 1 [over the suitable region where x goes from 0 to 2, causing K[sup]x[/sup] to trace out one full revolution]. In other words, the sum of K[sup]nx[/sup]/n over all nonzero integers n is constantly ln(K) * (1 - x).

Squaring the left hand side, it becomes a linear combination of K[sup]nx[/sup]s, each with coefficient equal to the sum of the reciprocal products of all pairs of nonzero integers adding up to n. In particular, when n = 0, this is the sum of the reciprocal products of all pairs of opposite nonzero integers; i.e., the sum of the reciprocal negated squares of the nonzero integers. In other words, the negation of our magic number.

Thus, if we could just find some linear operator which sends the function K[sup]nx[/sup] to 1 for n = 0 and to zero for other n, then our magic number would be equal to the absolute value of this operator on the function (ln(K) * (1 - x))[sup]2[/sup].

One particular such operator is the one which sends a function to its average value over the entire region of interest (from x = 0 to 2; for K[sup]nx[/sup] traces through n complete revolutions over this region, and thus, by symmetry, has average value zero for nonzero n). Our magic number is thus the absolute average value of (ln (K) * (1 - x))][sup]2[/sup] as x ranges from 0 to 2. [I.e., the variance of the value in radians of a random angle uniformly distributed through all 360 degrees]. That is, finishing off by “power rule” calculation, our magic number is |ln(K)|[sup]2[/sup]/3.

At this point, we note that, essentially by definition, |ln(K)| is π, and the argument is concluded.

57

A minor observation, perhaps of interest only to me: in fact, there is an equivalence here. The sum will cancel in the desired fashion if and only if K[sup]2c[/sup] - 1 = 0. However, that amounts to (K[sup]c[/sup] + 1)(K[sup]c[/sup] - 1) = 0, so, if K[sup]c[/sup] - 1 is to be invertible, the sum’s cancellation in the desired fashion is equivalent to K[sup]c[/sup] + 1 = 0.

58

(Thank you, Indistinguishable. I’ll try to find time to study your posts. I was already vaguely familiar with some of the relationship, and can even follow the proof of
http://mathworld.wolfram.com/EulerProduct.html
but it still all seems mysterious to me. :cool: )

Meanwhile, would this be a fun Golf problem? :
Given 4 (x,y) points, determine if they form a convex quadrilateral.

59

This is how I’d do it in Mathematica, but it is not very readable.


ConvexQuadQ[points_] := 
 And @@ (Union[{Sign[#[[1]][[1]]*#[[2]][[2]] - #[[1]][[2]]*#[[2]][[1]]], 
        Sign[#[[2]][[1]]*#[[3]][[2]] - #[[2]][[2]]*#[[3]][[1]]], 
        Sign[#[[3]][[1]]*#[[1]][[2]] - #[[3]][[2]]*#[[1]][[1]]]}] == 
        {-1, 1} & /@ Table[# - points[li] & /@ Delete[points, i], {i, 1, 4}])[/li]
60

Are the four points ordered, or should we consider all possible orderings of them?

In other words, would [ (0,0) , (0,1) , (1,0) , (1,1) ] be a convex quadrilateral?