noob python and C questions

[skeptical scientist]

I have a few noob python and C questions, and I expect to have more. Google is very helpful sometimes, but for questions I'm having trouble Googling, I was hoping to get some answers in this thread. Right now I just have one question, but I suspect I'll be having loads more, hence the generic thread title.

I've been working through the project Euler problems (currently in the low 60s), and I often find myself building some sort of data structure (e.g. a dictionary), and then iterating through it to find some particular result (the answer to the problem). If I need to iterate through the entire data structure, this is all well and good, but if I only care about iterating until I get the result (and then quit), I suspect some substantial gains in efficiency would be possible by only building the data structure as needed. (For example, I might only store a value for a given key in the dictionary when my search algorithm needs that particular value.)

Is there a graceful way to do this in Python? I can kludge it, but I suspect there's a much prettier way than my hack.

edit: added C questions

[thoughtfully]

Iterating through only part of a list is a classic use case for generators. Generators are a bit like functions that return a list of items, but one at a time, leaving you free to stop and go do something else if you like. You can save a lot of memory and usually speed, by not creating the whole list if you aren't going to the end anyway, or if the list is endless. There are a bunch of prime number generators out there, for instance.

A couple of math-nerd ones I made:
Returns the rows of Pascal's Triangle

Code: Select all

def pascal():
   r = 0
   l = [1]
   yield l

   while True:
      l  = [1] + [ l[i] + l[i+1] for i in range(r) ] + [1]
      r +=  1

      yield l


Returns the elements of the nth column of Pascal's Triangle. For n=2, this is the triangular numbers; for n=3, the tetrahedral numbers, and so forth.

Code: Select all

def triangular(n):
   if n == 0:
      while True:
         yield 1
   else:
      g = triangular(n-1)
      a = 0
      while True:
         a += g.next()
         yield a


Here's a good set of slides and sample code due to David Beazley, author of Python: The Essential Reference, the book I always have handy when coding. http://www.dabeaz.com/generators/

You can also Google "python generators" and get more resources than you could ever want.

You should probably also check out the Beginners Guide for Programmers, if you haven't. There's a similar guide for those new to coding.

Good luck with your adventures in Pythonland!

EDIT
A usage for you!

Code: Select all

# math.factorial() requires v2.6 or later
from math import log, factorial
from itertools import izip
import sys

# binomial coefficient
# also the nth row, kth column of pascal's triangle
# or n choose k
def bc(n,k):
   return factorial(n)/(factorial(n-k)*factorial(k))

def pascal():
   r = 0
   l = [1]
   yield l

   while True:
      l  = [1] + [ l[i] + l[i+1] for i in range(r) ] + [1]
      r +=  1

      yield l

def digits(n):
   return int(log(n)/log(10))+1

rows = int(sys.argv[1])

widths = [digits(bc(rows-1, i)) for i in range(rows)]

for i,row in izip(range(rows), pascal()):
   print ' '.join(['%%%dd'%(w,) for w in widths[:i+1]]) % tuple(row)


Code: Select all

0 % python pascal.py 10
1
1 1
1 2  1
1 3  3  1
1 4  6  4   1
1 5 10 10   5   1
1 6 15 20  15   6  1
1 7 21 35  35  21  7  1
1 8 28 56  70  56 28  8 1
1 9 36 84 126 126 84 36 9 1


[Ben-oni]

Actually, what you want are lazy data structures. This isn't directly available in Python, but it can be faked with lambdas, as in the following:

Code: Select all

def performLazyMap(f, xs):
   return map(lambda x: lambda: f(x), xs)

list = [0, 1, 2, 3, 4, 5]
lazyList = performLazyMap(lambda x: x*3, list)
print lazyList[3]() # outputs 3*3 = 9

Of course, the function still allocates the full array, which may not be desirable... if someone knows a more "Pythonesque" way to do the following, please enlighten me. I'm thinking of something along the C macro "#define CONS(x, y) cons(x, lambda: y)".

Code: Select all

# We need lazy lists, so we need to teach Python how to make those

# an empty list; also, list terminator
def null(b):
   return True if b==0 else null

# list creator
# cdr should always be wrapped in a lambda expression
def cons(car, cdr):
   forced = False
   def f(b):
      nonlocal forced
      nonlocal cdr
      if b < 0:
         return car
      elif b == 0:
         return False
      elif not forced:
         forced = True
         cdr = cdr()
      return cdr
   return f

def head(f):
   return f(-1)

def tail(f):
   return f(1)

def isNull(f):
   return f(0)

# We also need a convenient way to access members of a lazy list
def get(list, pos):
   while pos > 0:
      list = tail(list)
      pos = pos-1
   return head(list)

# And a few ways to work with a list would be nice
def lazyMap(f, xs):
   return cons(f(head(xs)), lambda: null if isNull(xs) else lazyMap(f, tail(xs)))

def zipWith(f, xs, ys):
   return cons(f(head(xs), head(ys)), lambda: null if isNull(xs) or isNull(ys) else zipWith(f, tail(xs), tail(ys)))

def iterate(f, x):
   return cons(x, lambda: iterate(f, f(x)))


naturals = iterate(lambda x: x+1, 0)

print(get(naturals, 1000))


As you see, it's possible to build infinite constructs, and do computations on them. You can even define them recursively:

Code: Select all

fibs = null
fibs = cons(1, lambda: cons(1, lambda: (zipWith(lambda x, y: x+y, fibs, tail(fibs)))))

print(get(fibs, 1000))


What this says is "The list of fibonacci numbers is equal to the list summed element-wise with the the list shifted over by one, given that the first two values are both 1." Functional declaration, linear time computation: the holy grail of programming.

EDIT: I suppose I'll throw in Pascal's, as well.

Code: Select all

def nextRow(row):
   def modZipWith(f, xs, ys, d):
      """The same as zipWith, except that we append a default value if any are missing"""
      x = d if isNull(xs) else head(xs)
      y = d if isNull(ys) else head(ys)
      return cons(f(x, y), lambda: null if isNull(xs) and isNull(ys) else (modZipWith(f, tail(xs), tail(ys), d)))
   return modZipWith(lambda x, y: x+y, cons(0, lambda: row), row, 0)

pascals = iterate(nextRow, cons(1, lambda: null))

def choose(n, m):
   return get(get(pascals, n), m)

print ("100 choose 50 =", choose(100, 50))


[skeptical scientist]

Yeah, generators look really cool. I haven't tried them for anything yet, but I suspect I will soon; just need the right application. I found a very pretty infinite prime generator using the sieve of Eratosthenes.

Generators seem great when you want to iterate over a list, and generate it as you go. But what if you want to use the list in more complicated ways? If you keep using the bit you've already generated, but occasionally you need to generate some new bit? What if instead of a list you want a mapping type, like a dictionary? What I have in mind is something like the following:

I have a function, f, which I need to compute for some large number of inputs, in such a way that I'm probably going to need to compute f on the same input some large number of times. Moreover, computing f is expensive. So to save time (at the expense of memory), I will use a dictionary to store computed pairs {n : f(n)}. Now, every time I want to compute f(n), I check whether n is a key in the dictionary. If it is, I return the corresponding value. Otherwise, I compute f, return f(n), and store the pair {n : f(n)} in my dictionary.

In certain applications, this could be much faster than precomputing the entire dictionary, and also much faster than computing f(n) every time you need it. I can do it with global variables, and I guess it's not too ugly that way, but I feel like there should be a nice local way of doing this.

[Ben-oni]

You're looking for memoization. In Python 3, it looks like this:

Code: Select all

def memoize(f):
   mem = {}
   def ret(x):
      nonlocal mem
      if not x in mem:
         mem[x] = f(x)
      return mem[x]
   return ret

# Let f be some big hairy function
memoizedF = memoize(f)


[Lothar]

Indeed, memoization seems to be what you want to do. Another way is for example:

Code: Select all

def fib(n, mem = {0:0, 1:1}):
    if n in mem:
        return mem[n]
    mem[n] = fib(n - 1) + fib(n - 2)
    return mem[n]

[Ben-oni]

Indeed, memoization seems to be what you want to do. Another way is for example:

Code: Select all

def fib(n, mem = {0:0, 1:1}):
    if n in mem:
        return mem[n]
    mem[n] = fib(n - 1) + fib(n - 2)
    return mem[n]

*shudder* Don't do that. It's not tail-recursive, and you will exceed your stack space.

[skeptical scientist]

*shudder* Don't do that. It's not tail-recursive, and you will exceed your stack space.

If it were tail-recursive, would that help? I didn't think Python optimized tail-recursions to not overflow the stack.

Here's a tail-recursive implementation of a gimped* Euclidean algorithm that overflows the stack (for me, running python 2.6.1 with mac defaults):

Code: Select all

def gcd(n,m):
    if m*n==0: return m+n
    if m>n:
        k=m
        m=n
        n=k
    return gcd(m,n-m)

print gcd(1000,1)


*Actually, this seems to be what Wikipedia calls the Euclidean algorithm. But in practice, you're much better off taking n mod m, rather than n-m.

[b.i.o]

I didn't think Python optimized tail-recursions to not overflow the stack.

It doesn't.

Another variation on what you were asking originally might be to implement lazy lists in Python. Creating a new class for them wouldn't be especially hard.

[EvanED]

I didn't think Python optimized tail-recursions to not overflow the stack.

It doesn't.

I was disappointed to see that this applies not just to CPython but to the JITting Pythons as well -- IronPython and PyPy. (I didn't try Jython and just assume it won't do it.)

I wonder if that's because the Python language makes guarantees (such as being able to access old stack frames in traceback or something) that makes tail call elimination hard.

[OmenPigeon]

That's basically it, except as far as I can tell it isn't guaranteed by the language, it's guaranteed by Guido. As I recall there was rather a kerfuffle in Python land around this point so I'm going to avoid wading into whether or not he's is in the right here, but the takeaway is that if you're want TCO in Python, your best bet is to outlive Guido.

[Ben-oni]

*shudder* Don't do that. It's not tail-recursive, and you will exceed your stack space.

If it were tail-recursive, would that help? I didn't think Python optimized tail-recursions to not overflow the stack.

Correct. Which is why I advocate lazy lists. Without tail-recursion, you will at some point require "trampolining", which is what my lazy-list indexer did (after a manner).

Here's a tail-recursive implementation of a gimped* Euclidean algorithm that overflows the stack (for me, running python 2.6.1 with mac defaults):

Code: Select all

def gcd(n,m):
    if m*n==0: return m+n
    if m>n:
        k=m
        m=n
        n=k
    return gcd(m,n-m)

print gcd(1000,1)


We can do this with trampolining:

Code: Select all

def t_gcd(n, m):
   if m*n==0: return lambda: m+n
   if m>n: return lambda: t_gcd(n, m-n)
   return lambda: t_gcd(m, n-m)

def trampoline(f):
   def g(*args):
      k = f(*args)
      while callable(k):
         k = k()
      return k
   return g

gcd = trampoline(t_gcd)

Once again, the implementation is lazy, but since the tail-recursion is being handled properly, calling "gcd(200, 2)" won't break the stack.

EDIT: A "better" way of doing this has occurred to me:

Code: Select all

def makeTailRecursive(funName):
   f = globals()[funName]
   def tF(*args): return lambda: f(*args)
   def g(*args):
      globals()[funName] = tF
      k = f(*args)
      while callable(k):
         k = k()
      globals()[funName] = g
      return k
   globals()[funName] = g

makeTailRecursive("gcd")
print(gcd(1000, 1))


Of course, it only works for top-level functions, but there may be a way to mess with the attribute stack... Someone care to give it a stab?

[EvanED]

That's basically it, except as far as I can tell it isn't guaranteed by the language, it's guaranteed by Guido.

Guido controls CPython, not PyPy etc.

If the language doesn't guarantee the absence of the tail call elimination (by means of guaranteeing a way to access backtraces that would make it "impossible" or some other way of making it "observable"), there's no reason other implementations couldn't do it; especially performance-minded implementations like PyPy.

Don't confuse implementation with language.

[OmenPigeon]

Guido controls CPython, not PyPy etc.

Sure, but it's not entirely that simple. As he mentioned in that blog post, tail calls are a feature, not an optimization. Adding them to a python implementation means that you can now write programs in one implementation of Python that won't execute correctly in others[1]. I'm not in the position of developing any Python fork, but if I were, I'd be really, really careful about letting my users write code that couldn't be used in CPython.

1: IronPython or Jython programs that rely on .NET or Java libraries also have this problem, but seems... less severe to me? Maybe because it's usually pretty obvious when you're using external libraries, but not when you're using fairly subtle language features?

[EvanED]

Guido controls CPython, not PyPy etc.

Sure, but it's not entirely that simple. As he mentioned in that blog post, tail calls are a feature, not an optimization.

This may be my fundamental disagreement.

I do think it's true that tail call optimization (TCO) is a feature in languages that rely on them heavily (functional languages that use recursion instead of loops), because otherwise recursion has very severe limitations.

However, that doesn't mean that can't ever be an optimization. For instance, GCC implements TCO for C, and this is an entirely reasonable decision. I doubt it guarantees it for C, and the C standard certainly doesn't. Remember, it's not just a matter of "now you have infinite stack frames and can execute stuff that wouldn't before"; it also improves memory and time use of your program.

If it were some complicated optimization I'd say it probably isn't worth it, but it doesn't seem like it is (modulo debug info); my only-slightly-naive position is that TCO seems like one of the easier optimizations to actually do.

Adding them to a python implementation means that you can now write programs in one implementation of Python that won't execute correctly in others[1]. I'm not in the position of developing any Python fork, but if I were, I'd be really, really careful about letting my users write code that couldn't be used in CPython.

TCO would make increase the difference sure, but it already is there. Different implementations already have different stack heights; this completes for me with PyPy 1.8 but not Python 2.7:

Spoiler:

Code: Select all

def countto(n, acc=0):
    if n == 0:
        return acc
    else:
        return countto(n-1, 1+acc)

print(countto(1000))

TCO just takes this same difference and makes it large.

And going the other direction, my impression is that there's a distressing amount of code that will work correctly in CPython but not in other implementations, because it implicitly relies on the reference-counting behavior of objects to do cleanup at the appropriate time.

This situation isn't ideal, but it's neither new nor made much worse by TCO.

[troyp]

Stackless Python does TCO. The only problem is that a handful of badly-behaved libraries don't work with it because they rely on the internals of interpreter. Unfortunately, one of them is PyQt (which is why I don't use Stackless, since I run Kubuntu and it would inevitably break something). Mind you, these days with virtualenv, it's easy to run multiple Pythons.

It really annoys me that Python lacks TCO, but Python (ie. Guido) is very anti-recursion in general^*, so you might as well just accept it. It's too painful trying to fight it.

Another thing (which is easily fixed, so it's not really annoying, just bizarre) is that Python seems to set the recursion limit insanely low. Tail-recursive code can usually be written naturally using iteration (hell, SICP even calls it iteration), but that's not always true for non-TC recursion, so I do sometimes use general recursion in Python. When I use a recursive algorithm on a problem of significant size, I always overrun the recursion limit and have to increase it. I don't quite understand why it's set so low in the first place.

* but what can you expect from someone who says that a generic for-loop is more readable (for experienced programmers, not noobs) than a reduce^**, and that the only use of reduce that's comprehensible is a trivial one with a commutative function.

** Thankfully, he supports list comprehensions (etc), which can often (though certainly not always) take the place of a fold^***. Which brings us back to Guido being an iteration man.
*** Reading this, it occurs to me this is completely untrue :-S I was talking about folds but thinking of higher order recursive functions in general. s/fold/higher-order list function/. I think I usually replace folds with either specialized folds (sum, max, etc)^**** or a for-loop. The for-loop often contains comprehensions to simplify it, though :-)
**** <pedantry> well, it's not really accurate to say I replace them, since I'd be using the specialized folds by preference in another language anyway </pedantry>

edit: re the original question, I agree that generators are going to be the main way to do lazy processing in Python. I also recommend generator expressions, which allow you to create generators inline in the same way you'd create a list with a list comprehension. For instance in Python 2.7

Code: Select all

for key in (key for key in mydict.keys()):
    # do stuff

lets you lazily iterate over the keys in a dict. In Python3, a lot of builtin functions, including dict.keys/values/items(), return generator-like objects, so you can just use

Code: Select all

for key in mydict.keys():
    # do stuff

Comprehensions and generator expressions are encouraged and well-optimized in Python and they can help make up for some of the things Python lacks, so they're a good thing to use. I find myself using them constantly.

[b.i.o]

List comprehensions are a replacement for map and filter, not fold. You can replace many uses of 'reduce' (which typically has the type 'a list -> ~f:('a -> 'a -> 'a) -> 'a) with a purpose-built function like 'sum' or 'max' (since those are both trivial to construct with reduce). Map and filter can only return lists. Reduce can only return a value of the type of thing stored in the list it's operating on.

Fold is different. You can't replace fold with anything but a loop (or recursion, obviously) in Python, because fold is far more powerful than any of the other standard higher-order functions. With filter or map you can only construct lists. With reduce you can only return a value of the original type. Fold's type is 'a list -> ~init:'b -> ~f:('b -> 'a -> 'b) -> 'b. 'b can be any type at all--fold is just as powerful as a 'for' loop with extra variables to act as accumulators.

(Python's reduce is actually a fold in disguise, I think, although I don't recall ever having used it. I separate the two because I think there's a clear distinction between the two and I think they're both useful functions in their own right, and it's much easier to be clear this way.)

[troyp]

Yeah, I was thinking of higher order functions in general. I corrected myself right afterwards but left the original text. You must have skipped over the correction.

edit:

(Python's reduce is actually a fold in disguise, I think, although I don't recall ever having used it. I separate the two because I think there's a clear distinction between the two and I think they're both useful functions in their own right, and it's much easier to be clear this way.)

reduce is just another name for foldl, except it has the init arg last (this is because init is optional, so Python's reduce actually doubles as fold and fold1).

I should add that reduce is equivalent to most versions of foldl. The scheme(r5rs/srfi1) version actually takes the arguments in reversed* order in its function argument. That is, reduce(f,xs,init) (and foldl f init xs in Haskell, etc) is equivalent to (foldl (lambda (a b) (f b a)) init xs) in Scheme.

* well, really it's everyone else who's using reversed order. I find this inconsistency really annoying, but the Scheme definition does make a lot of sense for any language that uses linked lists heavily. OTOH, they could have used the other definition and just added an "anticons" function for foldl-ing lists.

[skeptical scientist]

Noob C question:

I'm trying to write a program that does some simple command-line input and output. Sample code:

Code: Select all

int main() {
  int input=0;
  while(input<1 || input>10) {
    printf("Please enter a number between 1 and 10: ");
    scanf("%d", &input);   
  }
  printf("Your input was %d\n",input);
}

This works just as expected, as long as the user enters a number. If the number is in the specified range, the program spits back the number and terminates, otherwise it asks for a new number in the specified range. The problem is that if the user enters a non-integer, non-whitespace character, like a letter, the program goes into an infinite loop, because scanf keeps looking for a number, finding junk in the buffer, and returning without clearing the buffer. How can I clear the buffer at the beginning of the loop so this doesn't happen?

[EvanED]

See here. In short: read a line using fgets then parse it with sscanf. You can also clear until newline by using fgets and just ignoring the buffer.

BTW, sorry for the mini threadjack on TCO.

[skeptical scientist]

Thanks! That's perfect! It even has code which is identical to mine!

Random numbers are currently giving me a headache, because my google results are either giving me solutions that don't work on my system, or that give very poor results.

For example, what's a good way to generate a random number in a specified range? Some suggestions use expressions like

Code: Select all

rand()%N

to get a number in the range [0,N-1], but this clearly won't be uniform. (For example, if rand() returned numbers in the range [0,65535], and you used this method to get a random number in the range [0,40000], the numbers in 0-25534 would be each be twice as likely as numbers in 25535-40000.) Other suggestions use division, which seems more promising. I'm currently using

Code: Select all

rand()/(1+RAND_MAX/N)

which seems to give reasonable results in tests with ranges [0,1] and [0,19], but I'm not sure if the results are really uniform (to the extent that a pseudorandom generator can be called "uniform").

Also, how can I see my random number generator? The only suggestion I found online that works for my architecture was

Code: Select all

srand(time(NULL));

However, the results for this are laughably bad, since the first random number, at least, appears completely predictable.

[EvanED]

For example, what's a good way to generate a random number in a specified range? Some suggestions use expressions like

Code: Select all

rand()%N

to get a number in the range [0,N-1], but this clearly won't be uniform. (For example, if rand() returned numbers in the range [0,65535], and you used this method to get a random number in the range [0,40000], the numbers in 0-25534 would be each be twice as likely as numbers in 25535-40000.)

That's also not the greatest idea for more technical reasons; from what I've read, PRNGs tend to have higher entropy ("be more random") in the higher bits; taking the modulo leaves the low-order, lower entropy bits behind. (Or, at least it does if your range isn't relatively prime with RAND_MAX.)

What I do typically is

Code: Select all

rand() / (double) RAND_MAX * (high-low) + low

There's still a fundamental problem where it's mapping RAND_MAX integers onto a smaller range, so there will be numbers that are more probable than others, but they're not concentrated in one part of the range. If the (high-low) bit is small relative to RAND_MAX, the probability difference will be low as well.

If that still causes you consternation, you can find the largest multiple of (high-low) that is less than RAND_MAX, and if you draw a number that's larger, just discard it and try again. Then divide the result by (high-low).

Also, how can I see my random number generator? The only suggestion I found online that works for my architecture was

Code: Select all

srand(time(NULL));

However, the results for this are laughably bad, since the first random number, at least, appears completely predictable.

In what sense? It shouldn't be. If you run this a few times (at least a second apart), do you get different results?

Code: Select all

#include <stdlib.h>
#include <stdio.h>

int main() {
    srand(time(NULL));
    printf("%d\n", rand());
    return 0;
}

(Untested)

Edit: Oh, one common mistake is seeding too often. If you call srand() before each call to rand(), and you're using time(NULL), consecutive calls in the same process will often have the same seed and return the same value. Should still differ if you run it again though. This will almost always output two of the same value:

Code: Select all

#include <stdlib.h>
#include <stdio.h>

int main() {
    srand(time(NULL));
    printf("%d\n", rand());
    srand(time(NULL));
    printf("%d\n", rand());
    return 0;
}

If you need to get a seed more than once a second, you could also try seeding with some higher-frequency value, like the result from gettimeofday() on *nix, QueryPerformanceCounter() on Windows, or the rdtsc instruction on x86/x64.

[letterX]

However, the results for this are laughably bad, since the first random number, at least, appears completely predictable.

Also, if you happen to mean "predictable" in the security-sense of "if I knew what time you started your application, I could reconstruct all the random numbers you generate from that seed", then you're correct. If you need your random numbers to be unpredictable to outside observers, you should instead (at minimum) read your seed from /dev/random (at least if you're on *nix. Unsure about Windows). Also, at that point, you probably shouldn't use the standard rand() function, since some implementations of rand() can be pretty bad.

[skeptical scientist]

There's still a fundamental problem where it's mapping RAND_MAX integers onto a smaller range, so there will be numbers that are more probable than others, but they're not concentrated in one part of the range. If the (high-low) bit is small relative to RAND_MAX, the probability difference will be low as well.

I guess I can live with that, for most purposes. If I'm doing Monte-Carlo simulations, maybe your suggestion for better uniformity would be a good idea, but maybe I should just find a better PRG. If I'm doing crypto, I'll obviously have to find a real source of randomness.

Also, how can I see my random number generator? The only suggestion I found online that works for my architecture was

Code: Select all

srand(time(NULL));

However, the results for this are laughably bad, since the first random number, at least, appears completely predictable.

In what sense? It shouldn't be. If you run this a few times (at least a second apart), do you get different results?

I get different results for rand(), but they are similar. When mapping rand() to a D20 roll using the previously mentioned algorithm, I've gotten an 11 on the first roll every time I've tried it today. I would call that "laughably bad", and not in the sense of the previous poster.

[PM 2Ring]

Python stuff.

Generators seem great when you want to iterate over a list, and generate it as you go. But what if you want to use the list in more complicated ways? If you keep using the bit you've already generated, but occasionally you need to generate some new bit? What if instead of a list you want a mapping type, like a dictionary? What I have in mind is something like the following:

I have a function, f, which I need to compute for some large number of inputs, in such a way that I'm probably going to need to compute f on the same input some large number of times. Moreover, computing f is expensive. So to save time (at the expense of memory), I will use a dictionary to store computed pairs {n : f(n)}. Now, every time I want to compute f(n), I check whether n is a key in the dictionary. If it is, I return the corresponding value. Otherwise, I compute f, return f(n), and store the pair {n : f(n)} in my dictionary.

In certain applications, this could be much faster than precomputing the entire dictionary, and also much faster than computing f(n) every time you need it. I can do it with global variables, and I guess it's not too ugly that way, but I feel like there should be a nice local way of doing this.

List comprehensions are usually much faster than the equivalent for loop expression (at least in CPython), so a good compromise can be to pre-load the dictionary with some values that you expect you'll need.

You can use decorators to add memoization to a function. Of course, a decorator is just syntactic sugar, and you may prefer to create a memoized version of a function explicitly.

On the issue of local variables, remember that you can define functions inside functions in Python, and such functions have read access to the variables of the containing function. They can't modify those variables directly, but they can call their methods:

Code: Select all

>>> def a():
...  s = range(5)
...  def b():
...   s.append('z')
...  print s;b();print s
...
>>> a()
[0, 1, 2, 3, 4]
[0, 1, 2, 3, 4, 'z']


C stuff.
scanf() is evil. It's probably turned off more C newbies than any other feature of C. There are rarely any good reasons to use it.

To get decent random numbers in a given range out of rand() you'll need to use non-integer division. [ninja'd by EvanD]. I don't know why srand(time(NULL)) is misbehaving for you, however time() returns (time_t)(-1) if calendar time is not available.

The algorithm used by rand() isn't specified by the standard, and some old versions of rand() were pretty crummy. If you'd like to allow the user to specify random number seeds and you want to make your program portable, don't use the standard library's rand(), provide your own PRNG, preferably something modern like Mersenne Twister.

[skeptical scientist]

Python stuff.
On the issue of local variables, remember that you can define functions inside functions in Python, and such functions have read access to the variables of the containing function.

I didn't know that, actually. That's interesting.

C stuff.
scanf() is evil. It's probably turned off more C newbies than any other feature of C. There are rarely any good reasons to use it.

C generally seems much less friendly than Python. I coded a tic-tac-toe game in Python this morning and it took ~90 minutes. I've been working on a bagels game in C for more than 3 hours, and it's still only half-done. Granted, that's probably mostly because I've been solving project Euler problems in Python for the last two weeks, and I haven't really used C much in the last few years, but still...

The algorithm used by rand() isn't specified by the standard, and some old versions of rand() were pretty crummy. If you'd like to allow the user to specify random number seeds and you want to make your program portable, don't use the standard library's rand(), provide your own PRNG, preferably something modern like Mersenne Twister.

Right now I'm just making a stupid terminal game (bagels) so I don't need any strong randomness guarantees. I'd just like it if the bagels game didn't always choose the same "random" number!

[Jplus]

FYI: when you do need a better PRNG, it will be most convenient for you to look for a library. For example Boost.Random. It's C++, not C, but personally I think C++ would fit your needs better anyway. Otherwise, a PRNG library in pure C probably also exists.

Now that I mention it anyway, here's how your input/output problem would be solved in C++:

Code: Select all

#include <iostream>

int main() {
  using namespace std;
  int input=0;
  while(input<1 || input>10) {
    cout << "Please enter a number between 1 and 10: ";
    cin >> input;
  }
  cout << "Your input was " << input << endl;
  return 0;  // (added for correctness)
}

[EvanED]

FYI: when you do need a better PRNG, it will be most convenient for you to look for a library. For example Boost.Random. It's C++, not C, but personally I think C++ would fit your needs better anyway. Otherwise, a PRNG library in pure C probably also exists.

Now that I mention it anyway, here's how your input/output problem would be solved in C++:

Code: Select all

#include <iostream>

int main() {
  using namespace std;
  int input=0;
  while(input<1 || input>10) {
    cout << "Please enter a number between 1 and 10: ";
    cin >> input;
  }
  cout << "Your input was " << input << endl;
  return 0;  // (added for correctness)
}

FWIW, no it's not; that has the same problem as his original question. (Try entering something that isn't a number.)

Also, fun fact: in C++ you don't technically need the 'return 0;' at the end of 'main'; if execution falls off the end of main, it acts as it was there. ('main' is special that way; this doesn't apply to other functions, and doesn't apply to plain C.)

[Jplus]

Stuff[/code]

FWIW, no it's not; that has the same problem as his original question. (Try entering something that isn't a number.)

Oops, right. Here's the correct code:

Code: Select all

#include <iostream>

int main() {
    using namespace std;
    int input=0;
    while(input<1 || input>10) {
        cout << "Please enter a number between 1 and 10: ";
        cin >> input;
        if (cin.fail()) {
            cin.clear();
            cin.ignore(100, '\n');
            input = 0;
        }
    }
    cout << "Your input was " << input << endl;
}

Also, fun fact: in C++ you don't technically need the 'return 0;' at the end of 'main'; if execution falls off the end of main, it acts as it was there. ('main' is special that way; this doesn't apply to other functions, and doesn't apply to plain C.)

Thanks for pointing that out. So strictly speaking skeptical scientist was already coding in C++!

[Ben-oni]

Python stuff.
On the issue of local variables, remember that you can define functions inside functions in Python, and such functions have read access to the variables of the containing function.

I didn't know that, actually. That's interesting.

And in Python 3, write access as well.

Seriously, look at the code I provided for implementing lazy lists. It'll be worth the time spent trying to understand it. It uses: local definitions (nested def expressions), anonymous functions (lambda expressions), and scoped variables (the "nonlocal" keyword means you'll need to launch python3 to try it out, but you should have that available from Terminal). Once you understand those, see if you understand why accessing "fibs" to arbitrary depth works without overflowing the stack, or taking exponential time.

[skeptical scientist]

More noob Python Qs:

1. Is there any difference between x**0.5 and math.sqrt(x), and, if so, what is the difference?

2a. Can I can assume int(x**0.5) is the integer part of the square root of x so long as x < (sys.float_info.radix^{sys.float_info.mant_dig})²? If you check just the powers of 10, int(x**0.5) first fails to give the integer part of the square root of x for x=10³³, which is about what you would expect if this were true.

2b. If not, for what values of x can I safely assume that int(x**0.5) is the integer part of the square root of x? (Basically, if I care about precise integer arithmetic on long integers, when do I need to worry about writing my own integer square root function?)

3. Can anyone explain why the following script is causing my computer to run out of virtual memory:

Code: Select all

def int_sqrt(x): #returns the integer part of the square root of x, using a stupid algorithm. Only suitable near the limit of floating point accuracy.
    y=int(x**0.5)
    while y**2<x:
        y+=1
    while y**2>x:
        y-=1
    return y

for i in range(10000000):
    error=int(i**0.5)-int_sqrt(i)
    if error: print i,error

In order to execute this program, it seems to me, the only data that should be being stored is the current position in main, the current value of i, the current value of error, and, if int_sqrt is on the stack, the current position in int_sqrt, the current value of int_sqrt.x, and the current value of int_sqrt.y. It is true that the program calls int_sqrt some 10 million times, but never more than one instance is active at the same time. And yet, when I actually try to run the script, my system gets sluggish, and my available hard disk space drops by several hundred megabytes. I was expecting it to take ~10 minutes to execute; I was not expecting it to run my computer out of memory and force me to kill Python.

[letterX]

I don't know the answers to 1 and 2, and I'm actually really interested in knowing the answers myself. However for 3, have you tried changing

Code: Select all

for i in range(10000000):

to

Code: Select all

for i in xrange(10000000):

The thing about "range" is that it actually generates a list containing all numbers from 1 to n. As opposed to xrange, which just returns a generator which iterates through all numbers from 1 to n, without actually building the entire list.

[skeptical scientist]

For 3, have you tried changing

Code: Select all

for i in range(10000000):

to

Code: Select all

for i in xrange(10000000):

The thing about "range" is that it actually generates a list containing all numbers from 1 to n. As opposed to xrange, which just returns a generator which iterates through all numbers from 1 to n, without actually building the entire list.

Ah. Thanks, that sounds like a big improvement, since I don't want the whole list in memory at once.

[PM 2Ring]

1. Is there any difference between x**0.5 and math.sqrt(x), and, if so, what is the difference?

Well, math.sqrt(x) requires you to import the math module; apart from that AFAIK they may or not be identical, depending on the implementation. It's never a good idea to write code in any language that depends on quirks of the implementation.

If you want a Python integer sqrt that works reliably on large integers why not use Newton's method (aka Heron's method)?

[jareds]

These aren't really noob questions. In practice, a great number of experienced programmers seem to regard floating-point arithmetic as black magic. I don't endorse this attitude, but you may as well be aware of it. It looks like you already grasp the basics of floating-point representation.

1. Is there any difference between x**0.5 and math.sqrt(x), and, if so, what is the difference?

It is quite likely that x**0.5 calls the C library function pow(x,0.5), while math.sqrt(x) calls the C library function sqrt(x)¹. (If Python treats an exponent of 0.5 as a special case that also calls sqrt(x), then they are of course indistinguishable.) There are two potential differences between pow(x,0.5) and sqrt(x).

First, sqrt(x) is almost certain to be materially faster (which only matters if the square root computation is a bottleneck in the first place).

Second, sqrt(x) is required by C99 to be correctly rounded (as defined later), while pow(x,y) is not. As a practical matter, nowadays, I would say you have an extremely good chance that your C library has a correctly rounded sqrt, and at least a good chance that it has a correctly rounded pow; but this is basically an educated guess.

2a. Can I can assume int(x**0.5) is the integer part of the square root of x so long as x < (sys.float_info.radix^{sys.float_info.mant_dig})²? If you check just the powers of 10, int(x**0.5) first fails to give the integer part of the square root of x for x=10³³, which is about what you would expect if this were true.

First, there is the issue of whether x**0.5 or math.sqrt(x) is correctly rounded, the likelihood of which I just discussed. If it is correctly rounded, then it is the floating point number nearest to the mathematically exact result (if that is exactly halfway between two floating point numbers, then the correctly rounded result is the one of the two numbers whose mantissa ends in zero--this is called round-to-even). (Note that the mathematically exact result is completely unambiguous since x is simply a specific rational number.) IEEE actually defines other rounding modes, but when it comes to library functions, it's probably unwise to rely on any mode other than the one I just defined, also known as round-to-nearest. If it's not correctly rounded, then I would really expect it not to differ from the correctly rounded result by more than 1 or 2 ulps, which would then matter for large x, but if it's not correctly rounded then there's really nothing to formally rely on without a documented bound.

Second, we should consider what issues are possible if x**0.5 is correctly rounded. The potential problem is that y-0.5*ulp(y) <= square root of x < y for some integer y, in which case int(x**0.5) will be y when you want y-1. This is a legitimate potential problem. For large values in the range you gave, ulp(y) could be as large as 1.

Finally, I have been assuming all along that x is a floating-point number in the first place. If x is an integer that is not exactly representable as a floating-point number, you'll obviously need to that that into account, as x will first be converted to a floating-point number (by round-to-nearest unless you've changed the rounding mode).

I think at this point I have reduced this to a well-defined mathematical problem, which I am confident in leaving as an exercise for the reader.

2b. If not, for what values of x can I safely assume that int(x**0.5) is the integer part of the square root of x? (Basically, if I care about precise integer arithmetic on long integers, when do I need to worry about writing my own integer square root function?)

Personally, if x is an integer that is exactly representable as a floating-point number, then I would be perfectly confident doing int(math.sqrt(x)), and if I needed x to be larger than sys.float_info.radix^{sys.float_info.mant_dig}, I would probably go ahead and write an integer version without worrying about the exact bound, unless there were a reason for that bound to matter.

¹On x86, there is a hardware square root instruction, which the C library is basically certain to use. It's conceivable that Python also uses the hardware instruction directly without calling the C library, but this is then practically indistinguishable from calling the C library.

[skeptical scientist]

These aren't really noob questions.

I'm actually quite certain they are. I wrote helloworld.py on April 12, so I've known python for less than a month, which I'm pretty sure makes me a noob. And they are my questions. Hence, noob Python questions.

Anyways, thanks much for the responses to my questions. From your post, if x is an exact floating point representation of an integer (or an integer which has an exact floating point representation), then the only way int(sqrt(x)) can fail to be the integer part of the square root of x ([√x]) is when √x is slightly smaller than some integer n, but the closest floating-point representable number is larger than n, which means that int(sqrt(x)) will return n, but [√x]=n-1. As √x≤√(n²-1), the same issue will happen for n²-1, so to find the least point of failure, it suffices to search numbers which are one less than squares until an error is found, and then check between the failure point and the previous square to find the least failure point.

Here's a script that will find the least failure point:

Code: Select all

from math import sqrt

i=0
while 1:
    i+=1
    if int(sqrt(i**2-1))!=i-1: break
for j in xrange((i-1)**2,i**2):
    if int(sqrt(j))!=i-1: break
print "the least j for which int(sqrt(j)) is not the integer part of the square root of j is",j

On my system, it finds the failure point 4503599761588224, which is (2²⁶+1)²-1, after about 2 minutes. This makes sense, because √(n²-1) is approximately n-1/(2n), so √((2²⁶+1)²-1) should be approximately 2²⁶+1-1/2²⁷, which would require 54 bits of precision to accurately store (as 1.000 000 000 000 000 000 000 000 001 111 111 111 111 111 111 111 111 11₂*2²⁶). However, the square root of (2²⁶)²-1 would be approximately 2²⁶-1/2²⁷, which only requires 53 digits of precision to accurately store (as 1.111 111 111 111 111 111 111 111 111 111 111 111 111 111 111 111 111 1₂*2²⁵).

So it's not true that int(sqrt(x) returns [√x] for x<sys.float_info.radix^{sys.float_info.mant_dig}; however, it should be true in general that int(sqrt(x) returns [√x] for x<radix^mant_dig-1 when mant_dig is odd. I'm not quite sure where the point of failure should be when mant_dig is even, but surely not less than radix^mant_dig-2.

Since Python implements unbounded exact integer representations, it is possible that x may be an exact integer representation of a number which is not exactly representable as a float. In that case, int(sqrt(x)) might fail because x is first converted to a float, and so the result could be either larger or smaller than [√x], depending on what happens when x is converted to a float, and might be off by much more than just 1.

***

The main reason I'm asking this is I'm writing a lot of programs that need to find exact answers to mathematical questions (e.g. project Euler problems), I'm using int(sqrt(x)) a lot, sometimes on very large numbers, and I'm always a bit worried that I'm going to get the wrong answer because of inaccurate floating point arithmetic. I suspect I will end up needing to write my own exact int_sqrt() function at some point, I just want to know at what point I should bother.

[Meem1029]

Let's get one thing straight. No matter how long you have been programming, as soon as you are rightfully questioning implementation details about something like that, you are no longer asking noob questions.

[Jplus]

The main reason I'm asking this is I'm writing a lot of programs that need to find exact answers to mathematical questions (e.g. project Euler problems), I'm using int(sqrt(x)) a lot, sometimes on very large numbers, and I'm always a bit worried that I'm going to get the wrong answer because of inaccurate floating point arithmetic. I suspect I will end up needing to write my own exact int_sqrt() function at some point, I just want to know at what point I should bother.

I've solved quite a few Project Euler solutions, most of them in Python, and I've never run into problems with int(math.sqrt(x)). So I'd say: wait until you notice something isn't working.

By the way, I agree that your most recent questions are not really noob (well except for the range/xrange thing). You may be a Python noob, but you're obviously not a math noob.

[troyp]

I don't think it's inappropriate for skeptical to call himself a Python noob. I don't think he means to suggest that anyone who doesn't do error analyses of floating point functions is even more of a noob than him. It's just a matter of standards. It says "I don't think it's appropriate to call myself a competent Python programmer when I haven't even learned core language constructs yet." That's a good thing.

Of course his questions aren't all going to typical noob Python questions: he's a professional mathematician! Really, these are just the sort of questions you'd expect from a mathematician/computer scientist who was a Python noob.

[skeptical scientist]

Does anyone know of a nice python class/library for dealing with matrices?