Python, 27 bytes

lambda n:int((n*6)**.33359)

Try it online!

A direct formula with some curve-fitting, same as the original one found by Level River St.

The shifted equation i**3-i==n*6 is close to i**3==n*6 for large i. It solves to i=(n*6)**(1/3). Taking the floor the rounds down as needed, compensating for the off-by-one.

But, there are 6 inputs on boundaries where the error takes it below an integer it should be above. All of these can be fixed by slightly increasing the exponent without introducing further errors.

Python, 38 bytes

f=lambda n,i=1:i**3-i<n*6and-~f(n,i+1)

The formula n=i*(i+1)*(i+2)/6 for tetrahedral numbers can be more nicely written in i+1 as n*6=(i+1)**3-(i+1). So, we find the lowest i for which i**3-i<n*6. Each time we increment i starting from 1, the recursive calls adds 1 to the output. Starting from i=1 rather than i=0 compensates for the shift.

J, 12 bytes

2>.@-~3!inv]

There might be a golfier way to do this, but this is a lovely opportunity to use J's built-in function inversion.

Try it online!

How it works

2>.@-~3!inv]  Monadic verb. Argument: n

           ]  Right argument; yield n.
      3       Yield 3.
       !inv   Apply the inverse of the ! verb to n and 3. This yields a real number.
              x!y computes Π(y)/(Π(y-x)Π(x)), where Π is the extnsion of the 
              factorial function to the real numbers. When x and y are non-negative
              integers, this equals yCx, the x-combinations of a set of order y.
 >.@-~        Combine the ceil verb (>.) atop (@) the subtraction verb (-) with
              swapped arguments (~).
2             Call it the combined verbs on the previous result and 2.

Python, 22 bytes

lambda n:n**.3335//.55

Heavily inspired by @xnor's Python answer.

Try it online!

Jelly, 7 bytes

R‘c3<¹S

Try it online!

How it works

R‘c3<¹S  Main link. Argument: n

R        Range; yield [1, ..., n].
 ‘       Increment; yield [2, ..., n+1].
  c3     Combinations; yield [C(2,3), ..., C(n+1,3)].
    <¹   Yield [C(2,3) < n, ..., C(n+1,3) < n].
      S  Sum; count the non-negative values of k for which C(k+2,3) < n.

JavaScript (ES6), 33 bytes

n=>(F=k=>k<n?F(k+3*k/i++):i)(i=1)

Based on the recursive formula:

a(1) = 1
a(i) = (i + 3) * a(i - 1) / i

The second expression can also be written as ...

a(i) = a(i - 1) + 3 * a(i - 1) / i

... which is the one that we are using here.

a(i - 1) is actually stored in the k variable and passed to the next iteration until k >= n.

Test cases

let f =

n=>(F=k=>k<n?F(k+3*k/i++):i)(i=1)

console.log(f(1));   // ->  1
console.log(f(5));   // ->  3
console.log(f(75));  // ->  7
console.log(f(100)); // ->  8
console.log(f(220)); // ->  10
console.log(f(221)); // ->  11
console.log(f(364)); // ->  12

JavaScript, 36 33 bytes

-3 bytes thanks to Luke (making the function curried)

n=>f=i=>n<=i/6*-~i*(i+2)?i:f(-~i)

This is an unnamed lambda function which can be assigned to func and called with func(220)(), as described in this meta post. My original, non-curried function looks like this:

f=(n,i)=>n<=-~i*i/6*(i+2)?i:f(n,-~i)

This answer uses the fact that xth tetrahedral number can be found with the following function:

The submission works by recursively increasing i, and finding tetrahedral(i), until it's larger than or equal to n (the number of presents given).

When called with one argument as expected, i = undefined, and therefore is not larger than n. This means f(n,-~i) is executed, and -~undefined evaluates to 1, which sets off the recursion.

Test Snippet:

func = n=>f=i=>n<=i/6*-~i*(i+2)?i:f(-~i)

var tests = [1, 5, 75, 100, 220, 221, 364];
tests.forEach(n => console.log(n + ' => ' + func(n)()));

Ruby, 26 bytes

Edit: alternate version, still 26 bytes

->n{(n**0.3333*1.82).to_i}

Original version

->n{((n*6)**0.33355).to_i}

Uses the fact that T(x) = x(x+1)(x+2)/6 = ((x+1)**3-(x+1))/6 which is very close to (x+1)**3/6.

The function simply multiplies by 6, finds a slightly tweaked version of the cube root (yes 5 decimal places are required) and returns the result truncated to an integer.

Test program and output

f=->n{((n*6)**0.33355).to_i}
[1,4,10,20,35,56,84,120,165,220,286,364].map{|i|p [i,f[i],f[i+1]]}

[1, 1, 2]
[4, 2, 3]
[10, 3, 4]
[20, 4, 5]
[35, 5, 6]
[56, 6, 7]
[84, 7, 8]
[120, 8, 9]
[165, 9, 10]
[220, 10, 11]
[286, 11, 12]
[364, 12, 13]

Jelly, 7 6 bytes

-1 byte thanks to Dennis (use vectorised minimum, «, and first index, i)

R+\⁺«i

TryItOnline

How?

Not all that efficient - calculates the 1st through to nth tetrahedral numbers in order in a list and returns the 1-based index of the first that is equal to or greater.

R+\⁺«i - main link: n
R      - range                          [1,2,3,4,...,n]
 +\    - cumulative reduce by addition  [1,3,6,10,...,sum([1,2,3,4,...n])] i.e. triangle numbers
   ⁺   - duplicate previous link - another cumulative reduce by addition
                                        [1,4,10,20,...,nth tetrahedral]
    «  - min(that, n)                   [1,4,10,20,...,n,n,n]
     i - first index of n (e.g. if n=12:[1,4,10,12,12,12,12,12,12,12,12,12] -> 4)

Previous 7 byters using lowered range [0,1,2,3,...,n-1] and counting tetrahedrals less than n:
Ḷ+\⁺<µS,
Ḷ+\⁺<ḅ1,
Ḷ+\⁺<ċ1, and
Ḷ+\⁺<¹S

MATL, 12 11 bytes

`G@H+IXn>}@

Try it online!

Explanation

`       % Do...while
  G     %   Push input, n
  @     %   Push iteration index (1-based), say m
  H     %   Push 2
  +     %   Add
  I     %   Push 3
  Xn    %   Binomial coefficient with inputs m+2, 3
  >     %   Is n greater than the binomial coefficient? If so: next iteration
}       %   Finally (execute after last iteration, before exiting the loop)
  @     %   Push last iteration index. This is the desired result
        % End (implicit)
        % Display (implicit)

05AB1E, 10 bytes

XµNÌ3c¹‹_½

Try it online!

Explanation

Xµ          # loop until counter equals 1
  NÌ3c      # binomial_coefficient(iteration_index+2,3)
      ¹     # push input
       ‹_   # not less than
         ½  # if true, increment counter
            # output last iteration index

Pyke, 11 bytes

#3RL+B6f)lt

Try it here!

#3RL+B6f)   -  while rtn <= input as i:
 3RL+       -     i+j for j in range(3)
     B      -    product(^)
      6f    -   ^/6
         lt - len(^)-1

Mathematica, 26 bytes

(0//.i_/;i+6#>i^3:>i+1)-1&

Unnamed function taking a nonnegative integer argument and returning a nonnegative integer (yeah, it works for day 0 too). We want to find the smallest integer i for which the input # is at most i(i+1)(i+2)/6, which is the formula for the number of gifts given on the first i days. Through mild algebraic trickery, the inequality # ≤ i(i+1)(i+2)/6 is equivalent to (i+1) + 6# ≤ (i+1)^3. So the structure 0//.i_/;i+6#>i^3:>i+1 starts with a 0 and keeps adding 1 as long as the test i+6#>i^3 is satisfied; then (...)-1& subtracts 1 at the end (rather than spend bytes with parentheses inside the inequality).

If we let the 12 Days of Christmas continue, we can handle about 65536 days before the built-in recursion limit for //. halts the process ... that's about 4.7 * 10^13 days, or about ten times the age of the universe thus far....

J, 9 bytes

I.~3!2+i.

Try it online!

This is more inefficient than using the inverse of factorial but happens to be shorter.

For example, if the input integer is n = 5, make the range [2, n+1].

2 3 4 5 6 choose 3
0 1 4 10 20

These are the first 5 tetrahedral numbers. The next step is to determine which interval (day) n belongs to. There are n+1 = 6 intervals.

0 (-∞, 0]
1 (0, 1]
2 (1, 4]
3 (4, 10]
4 (10, 20]
5 (20, ∞)

Then n = 5 belongs to interval 3 which is (4, 10] and the result is 3.

Explanation

I.~3!2+i.  Input: integer n
       i.  Range [0, n)
     2+    Add 2 to each
   3!      Combinations nCr with r = 3
I.~        Interval index of n

Python, 43 bytes

f=lambda n,i=0:n*6>-~i*i*(i+2)and-~f(n,i+1)

Saved 5 bytes thanks to @FlipTack and another 3 thanks to @xnor!

SmileBASIC, 43 bytes

INPUT X
WHILE X>P
I=I+1
R=R+I
P=P+R
WEND?I

I is the day, R is the ith triangular number, and P is the ith tetrahedral number (number of presents).

I think a similar answer in another language, perhaps: x=>{while(x>p)p+=r+=++i;return i} could be pretty good.

Python 3, 48 46 bytes

f=lambda x,i=1:f(x,i+1)if(i+3)*i+2<x/i*6else i

Haskell, 21 bytes

floor.(/0.82).(**0.4)

Only works on the exact range of 1 to 12 days. You can test with

λ map(floor.(/0.82).(**0.4))$map(\n->sum$(reverse>>=zipWith(*))[1..n])[1..20]

[1,2,3,4,5,6,7,8,9,10,11,12,14,15,16,17,19,20,21,22]

Mathematica, 31 25 bytes

Floor@Root[x^3-x-6#+6,1]&

Japt, 12 bytes

1n@*6§X³-X}a

13-byte solution, inspired by @xnor's answer:

p.3335 /.55 f

A few more solutions @ETHproductions and I played around with

J+@*6§X³-X}a 
@*6§X³-X}a -1
@§X/6*°X*°X}a 
_³-V /6¨U}a -1
§(°V nV³ /6?´V:ß
§(°VV³-V /6?´V:ß

Test it here.

Batch, 69 bytes

@set/an=d=t=0
:l
@set/at+=d+=n+=1
@if %t% lss %1 goto l
@echo %n%

Manually calculates tetrahedronal numbers.

Pyth 11 bytes

/^Q.3335.55

Try it online!

pretty much just translated Dennis' answer into Pyth

R, 19 characters

floor((n*6)^.33359)

based on xnor's answer in Python.

QBIC, 19 bytes

This steals @xnor 's formula:

:?int((a*6)^.33359)

I tried dialing down the resolution to .3336 to save a byte, but that fails on the final testcase.